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Abstract 



This appendix for the author's forthcoming thesis, "On-Line Solu- 
tion of Elliptic Boundary- Value Problems," is a user's guide for EPS. 
EPS solves two-dimensional boundary-value problems for elliptic systems 
of second-order partial differential equations. It also has general- 
purpose capabilities which permit the on-line definition and execution 
of arbitrary numerical procedures. 

The guide is concerned primarily with the use of EPS for solving 
elliptic boundary- value problems. Linear problems of this type that 
have no complications such as free surfaces or undefined parameters can 
be solved on a one-pass basis. Nonlinearities and other complications 
can be accommodated by iteration. Solutions are obtained by a finite- 
difference method which permits the use of irregular lattices, hence the 
crowding of nodes in sensitive regions. 

EPS operates on the IBM 7094 computer of the M.I.T. Compatible 
Time-Sharing System (CTSS) , and exploits to an unusual degree the poten- 
tial for interactive problem solving that CTSS affords. Input commands 
resemble statements in various algebraic compiler languages, and can be 
combined and abbreviated by means of macros. Improper input and other 
error conditions are handled so as to minimize user inconvenience. 
Common syntax errors, for example, are corrected automatically by the 
machine. Output is available in either numerical or graphical form. 
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Preface 

This user's guide for EPS will be included as an appendix in the 
thesis "On-Line Solution of Elliptic Boundary-Value Problems" [1]. How- 
ever, for the convenience of current and potential EPS users, it is be- 
ing made available in advance of the entire thesis as a separate, self- 
contained document. As this document necessarily specifies the capabil- 
ities and important features of EPS, it may be of interest not only to 
actual users, but also to others with interests in such areas as man- 
machine communication, problem-oriented language design, and numerical 
analysis . 

The guide has sixteen major sections, the last being a summary of 
the preceeding fifteen. For purposes of orientation, it may be helpful 
to read the summary first (p a g es 168-177) . 

Aside from the summary, many users will find that only the first 
seven sections are essential. Sections 1 and 2 discuss program access 
and basic operating conventions, while Sections 3 through 7 discuss the 
steps required for solving two-dimensional boundary-value problems for 
linear elliptic systems of second-order partial differential equations. 
Users with problems of this class that are initially well-defined, i.e., 
that have no free surfaces or parameters needing optimization, should m 
principle be able to ignore the rest of the guide, although operation of 
the system may be simplified by using features discussed later on. 



On the other hand, users with problems that have nonlinearities or 
other complications, as well as those who wish to investigate various 
special extensions, e.g., the generation of flow nets or the use of EPS 
as a general-purpose language, will probably find the rest of the guide 
essential, particularly Sections 8, 9, and 10. 

All users may need to refer occasionally to Sections 12 through 15 
in order to overcome various possible difficulties. 
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Introduction 

When the development of EPS first began in the fall of 1964, it 
was intended that this system be a fully automatic, numerical assistant 
which would accept "textbook problem descriptions" of a particular class 
and produce reliable answers without demanding an intimate knowledge of 
computer methods from its users. Specifically, it was intended that EPS 
solve two-dimensional boundary-value problems for linear elliptic sys- 
tems of second-order partial differential equations, and it was hoped 
that users would be able to obtain solutions to particular problems by 
supplying only essential information describing governing equations, 
boundary geometry, and boundary conditions. 

In a sense, EPS has both overshot and fallen short of its original 
goals. It has overshot these goals insofar as it is able today to treat 
a class of problems which is much broader than, although not so rigor- 
ously definable as, the class for which it was initially intended. On 
the other hand, EPS is not the "fully automatic, numerical assistant" 
once envisioned. With the benefit of hindsight, it is possible to say 
that both of these results were inevitable. The class of problems which 
EPS was meant to solve is itself broad enough to indicate the use of 



EPS is an acronym for "Equilibrium Problem Solver," an expression 
suggested by the fact that in mathematical physics boundary-value prob- 
lems generally arise when continuum models are employed in the analysis 
of physical systems at equilibrium. 
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numerical methods which, for maximum flexibility, are best implemented 
in a semi-automatic form. Thus, in particular, it is necessary in EPS 
for the user to define a discrete model for his problem, in the form of 
a finite-difference lattice, in addition to the various parameters which 
are mathematically essential. At the same time, however, the many vari- 
ables left open by the class of problems for which EPS was intended ne- 
cessitated certain general-purpose, algebraic capabilities which make 
possible many important extensions. For example, it is possible to have 
EPS compute revised problem or lattice parameters using results from 
previous solutions. Consequently, iterative procedures can be effected 
which permit the treatment of boundary-value problems with nonlineari- 
ties, free surfaces, or undefined parameters. Further, the general-pur- 
pose features of the system can be exploited in an independent manner, 
and enable the user to specify and have executed algorithms which have 
nothing to do with the solution of boundary-value problems. 

While in principle the various capabilities of EPS derive solely 
from specific computational features, in practice the utility and use- 
ability of the system stem as well from the fact that it was designed 
for interactive operation in an on-line computing environment. This 
arrangement is especially convenient, perhaps essential, in sophisti- 
cated applications, where the user may need to make numerous decisions 
in order to obtain solutions. Moreover, as the beginning user will 
quickly appreciate, the task of learning proper operating procedures 
is greatly facilitated when the results of a particular action can be 
determined immediately. In the design of EPS, considerable effort has 
been made to insure that input errors are detected quickly and reported 
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in an instructive manner, and, in general, there is no danger of des- 
troying the system or losing past results because of erroneous input. 
Therefore, the user should not only feel free to experiment, but should 
realize that experimentation may well be the most expedient way of deter- 
ning whether a specific input form or solution strategy is valid. 



mi 
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1. Access 

EPS was implemented for operation on the IBM 7094 computer of 
the M.l.T. Compatible Time-Sharing System (CTSS) . It is assumed that 
the reader has access to CTSS and, further, that he is familiar with 
basic CTSS operating procedures and nomenclature. (The beginning user 

may find the "Time-Sharing Primer" in [2], Section AA.2, to be help- 

. f 
ful in this regard.) Meeting these requirments, he may log m at 

a computer console and proceed as indicated below. 

There are currently two versions of EPS. One is maintained in 

a disk file called "eps saved" and is intended for users who have no 

desire, either present or anticipated, for graphical output. The other 

is maintained in a file called "eps-p saved" and is intended for users 

who do wish graphical output and who therefore have either logged in 

§ 
at an ESL graphic display console or anticipate continuing their com- 
putations from such a console at a later time. To establish a link 

11* 
to either file, type 



link name saved mlAl6 cmfl04 



where name is the appropriate first file name, i.e., "eps" or "eps-p 



See [2], sees. AA.2 and AH. 1.01. 

See [2], sees. AC.O and AH. 2.06, and Thornhill, et al [3]. 

EPS is available through the CTSS public file directory. Cf. 

[2], sees. AD. 4 and AH. 3. 05. 
* 
This guide uses lower-case letters for computer input and upper- 
case letters for computer output, as is the normal convention on most 
CTSS consoles. Characters intended literally are always shown in sans- 
serif type, while variable information is shown in italics or other 
easily distinguished type styles. 



The computer will respond with lines of the form 

W T 

R A5 n 1 +A.T 2 

where T, AT-^ and A2 7 2 represent timing information, as explained in [2] 

t§ 
After this, type 



resume name^ 



to initiate a run. The computer will respond by printing 



W T 

THIS VERSION OF EPS WAS LAST REVISED ON DATE. 

PROCEED: 



With this, submission of input can begin. 



Cf . [21 , sec. AH. 7. 03. 



If the user is stationed at and planning to use an ESL display, 
he may need to load the PDP-7 display monitor before resuming EPS. See 
posted instructions or ESL staff for current procedures. 



2. Basic operating conventions 

Like CTSS, EPS performs the various operations of which it is 
capable in response to explicit commands from the user. Commands can 
be supplied directly from the console or they can be supplied from a 
disk file prepared in advance. However, it is recommended that only 
console input be used initially, if not in general, as it lends itself 
to faster understanding and poses fewer complications when mistakes are 
made. Therefore, most of this guide, including the present section, 
assumes that input is being typed in directly. Disk input, which in 
any event is basically simulated console input, is discussed further in 
Section 11, as is a mechanism for obtaining disk output. 

In contrast to CTSS, the basic unit of input to EPS is not neces- 
sarily a single line of type representing a single command. Instead, 
any number of lines may be typed, and these may represent any number 
of commands. (No correspondence between physical lines and individual 
commands is required: one command may follow another on the same line, 
or a command may be continued from one line to the next, without special 
punctuation or continuation marks.) Thus the general procedure for 
using EPS is as follows: After receiving a "PROCEED:", the user types 
one or more commands on as many lines as are necessary and convenient, 
the maximum line length being 83 characters. After typing the last 
command and before giving a final carriage return, he types a "$" to 
indicate that processing can commence. The computer executes the com- 
mands in order of their specification, then prints another "PROCEED:" 
to indicate that more input can be submitted. This procedure is illus- 



trated by the following sequence, which happens to show the definition 
of a finite-difference lattice for a semi-annular domain: 



PROCEED: 

append 0,0, 8,0, 8,20, 0,20 to zeta 

close zeta$ 

DEFINITION OF NEW CURVE 'ZETA' HAS BEEN COMPLETED. 

'ZETA' HAS BEEN CLOSED. 

PROCEED: 

define x-radius*cos(theta) , y=radius*sf n(theta) , 

radius=ri+(ro-ri)*l/8, theta-pi*j/20 
set ri-2.2, ro-4.7, pl-4*atan(l) , \-k, j-0 
print x,y $ 

3.4500000E+00 0.00000000000 

PROCEED: 

eta. 



It is absolutely essential that the user type a "$" and a_ carriage re- 
turn on completing each unit of input; otherwise, the computer will wait 
for more input and accomplish nothing. Eventually, of course, the user 
will wish to stop interacting with EPS and go on to other matters. To 
do this, he can type 

quit$ 

whereupon control will be returned to the CTSS supervisor. For further 
details, see Section 13. 

As the examples above indicate, commands for EPS are basically 
like statements in various algebraic compiler languages. However, cer- 
tain common restrictions have been eliminated in the interest of more 
natural and flexible communication. The following features should be 
noted, in particular: 
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(1) An identifier, e.g., the name of an algebraic parameter, may be 

any string of up to 83 letters and digits that begins with a let- 

f 
ter and is not a reserved word. 



(2) A number may be written in integer, decimal, or scientific (E-) 
format, regardless of the type of data it represents. Thus the 
following are equivalent: 30, 30.0, 3.0e+01, 3el, .3e2, 300e-l. 

(3) Command arguments that represent numerical values may be written 
as algebraic expressions in terms of numbers, parameter names, 
built-in functions, and standard algebraic, relational, and Bool- 
ean operators. Further, because Boolean values are represented 
internally as numerical l's and 0's, and because no parameter mode 
distinctions are observed, it is permissible to mix operations of 
different types at will. For example, the expression 

(s grt 5 and s les 7-9)"t 

has the value "l»t -> t" when "s" lies between 5 and 7.9, and the 
value "0 ; ' : t -> 0" otherwise. 

(4) Explicit storage declarations for subscripted parameters, etc., 
are unnecessary, as storage is allocated and reallocated automa- 
tically as problem information is introduced and deleted. 



4- 

Reserved words are the command words, multi literal operators, and 
built-in function names shown in sans-serif type in Tables 1-3 of sec. 
16, plus "if", "insert", and "bed". The symbol "bed" is used to intro- 
duce literals, which, though intended primarily for specifying output 
formats, may also be used as identifiers. See sees. 7.1 and 15.1. 

See sec. 15.1 and Tables 2 and 3, sec. 16. 
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Other conventions are more or less standard. For example, blanks are 
optional between any input element and any printed delimiter , namely, 
any one of the following: 

+ -/*() = : , $ 

On the other hand, at least one blank must be introduced between adja- 
cent words and numbers. Tabs and carriage returns are considered to be 
logically equivalent to blanks; therefore, it is not possible to split 
up a word or number, typing part on one line and the rest on the next. 
Typing errors in individual lines of EPS input can be corrected 
by the same methods used to correct command lines to CTSS, namely, by 
typing a question mark (?) or at-sign (@) to cause all preceeding char- 
acters in a given line to be ignored, or by typing one or more quota- 
tion marks (") or number signs (#) to cause one or more immediately 
preceeding characters to be ignored. Further, the processing of input, 
once it has been typed, can be curtailed, if not prevented entirely, 
by use of the program interrupt button (see Section 13), while errone- 
ous data can always be corrected by typing new commands which in effect 
override old ones (see Section 8). In any event, as indicated earlier, 
errors which are not detected by the user are generally detected and 
agreeably reported by the computer; in fact, many common syntax errors 
are not only detected but also corrected automatically (see Section 14). 



Beginning users should beware that, because EPS uses the CTSS 
6-bit character set, certain symbols available on some consoles, e.g., 
the [ ] % C ! etc., are ignored. Thus typing sin(t%") is equivalent 
to typing sin() , not_ sin(t) . Cf . [2] , sec. AC . 2 . 



3. Lattice specification 

Input required by EPS for the solution of a two-dimensional bound- 
ary-value problem can be divided into three catagories: 

(1) information describing the finite-difference lattice by 
which the domain of interest is to be modeled; 

(2) information describing the governing partial differential 
equation (s) for the domain; and 

(3) information describing the auxiliary conditions to be 
imposed along the contour (s) which bound the domain. 



on 
com- 



Since EPS does not actually assemble such information until a soluti 
is explicitly requested, the order in which it is presented to the 
puter is virtually arbitrary. However, for the sake of concreteness , 
it is assumed here that the user will specify his input in the order 
indicated above. Therefore we first consider the selection and speci- 
fication of a finite-difference lattice. 

3.1. Admissible configurations. In the analysis of two-dimen- 
sional boundary-value problems by conventional finite-difference methods 
it is common to employ lattices based on orthogonal grids of equally 
spaced, straight lines. In some cases the analysis is carried out by 
writing and solving difference equations for those points which lie 
within or on the boundary of the domain of interest and which repre- 
sent intersections of the grid lines with one another or with the bound- 
ary (see Figure la). The admission of boundary points which do not 
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Conventional lattices based on regular, orthogonal 
grids. Labeling conventions employed with EPS 
prohibit lattice (a)^ but admit lattice (b). 



coincide with the nodes of the grid, i.e., the points at which the grid 
lines intersect one another, permits accurate modeling of boundary geom- 
etry, but introduces the possibility of irregularities in lattice topol- 
ogy which complicate automatic analysis. For this reason, many computer 
programs require that domain boundaries be represented (in general, 
approximated) by polygons whose sides join neighboring or diagonally 
opposed nodes (see Figure lb). In fact, a requirement of this type is 
imposed by EPS, so that the lattice of Figure lb is admissible to EPS 
while that of Figure la is not. Fortunately, however, Figure lb is not 
representative of the entire class of lattices which EPS will accept, 
for EPS does not demand that grid lines be orthogonal, equally spaced, 
or straight. (It assumes that nodes are connected with straight lines, 
but allows these lines to vary in length and direction from node to 
node.) Consequently, the domain of Figure 1 could also be modeled with 
any of the lattices shown in Figure 2. 
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Geometrically irregular lattices, all admissible 
to EPS, with possible grid coordinate labeling s. 



It is possible to label all lattice points in Figure lb, and 
also all points in the lattices of Figure 2, in a discrete manner, 
assigning each node a unique pair of integers, or grid coordinates , 
such that the immediate neighbors, listed in clockwise order, of any 
node (i,j) not on a lattice boundary are (i+l,j) , (i,j-l), (i-l,j), 
and (i,j+l). (Note that such a labeling is impossible for Figure la 
because of the special points along the curved portion of the boundary.) 
With this labeling convention, grid lines can be defined formally as 
paths through a lattice connecting nodes of constant i or constant j. 
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Clearly, if such a labeling is possible at all, then it can be carried 
out in any number of ways. However, once an -£- or j -direction and the 
grid coordinates for any single node have been selected, the coordi- 
nates of all other nodes are fixed, and can be determined by counting. 
It can be stated categorically that any finite lattice which can 
be labeled with discrete grid coordinates in the manner described above 
is, in principle, admissible to EPS. Conversely, lattices that cannot 
be so labeled are not admissible. In order for an acceptible labeling 
to be possible, a lattice must satisfy the following requirements: 

(1) All lattice points must represent the intersections, or nodes, 
of two families of grid lines (lines of constant i and lines 
of constant j), arranged in such a way that each node not on 
the lattice boundary is connected to four and only four neigh- 
bors. This requirement specifically prohibits grid lines 
which terminate at interior nodes, a possibility sometimes 
admitted to allow for abrupt changes in lattice spacing. 

(2) Neighboring grid lines of the same family, i.e., lines of 
constant i (or j) with i- (or j-) values differing by ±1, 
may not coincide with or cross over one another. However, 
in principle at least, several lines of the same family may 
converge to a single point (and then diverge away from that 
point, if the point is not on the lattice boundary), produc- 
ing a set of superimposed nodes. (See Section 3.2 for prac- 
tical recommendations.) 



•10- 



[3) The lattice boundary must consist of one or more closed poly- 
gons, the sides of which either coincide with grid lines or 
form grid diagonals , i.e., lines connecting nodes whose grid 
coordinates satisfy a relation of the form 



integer 
constant 



Note that in general any admissible lattice will appear as a net of 
quadrilateral and triangular elements, the latter arising (if at all) 
either along boundary segments which represent grid diagonals or in 
places where quadrilateral element degenerate into triangles because 
of the use of superimposed nodes. 

3.2 Practical limitations and recommendations. In principle, 
the number of admissible lattices for any given problem is unlimited, 
not only because the number of grid lines and nodes is variable, but 
also because there are frequently several choices as to which basic 
form a lattice should assume. In practice, of course, there are vari- 
ous limiting considerations, the most important of which are discussed 
below. 

(1) For the following reasons, it is generally impractical, if 
not impossible, to employ lattices with more than a few hun- 
dred nodes : 

(a) The generation of finite-difference equations for a 

system of n second-order partial differential equations 
and a lattice with m nodes requires about n z m/10 sec- 
onds of IBM 7094 computer time. Time requirements 
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for the solution of these equations and for the gener- 
ation of output vary considerably, but taken together 
are typically of the same order. 

(b) The number of words of core storage W required to treat 
a particular problem varies with 

m total number of nodes 

b number of boundary nodes 

n number of partial differential equations 

nfa number of dependent variables in fcth equation 

d number of distinct sets of difference equa- 
tion coefficients (see below) 

p number of storage words needed to represent 
algebraic parameters and macros used in prob- 
lem definition (typically 500 to 2000 words) 



and is given approximately by 

W * rnn + 8b + 



n 
l + 9ln k d + p 
k=l > 



Clearly W must not exceed the amount of storage actu- 
ally available, namely, 20,500 words with "eps saved" 
and 18,600 words with "eps-p saved". For regular lat- 
tices and differential systems with constant coeffi- 
cients, it is likely that only a few distinct sets of 
difference equation coefficients will be stored {d = 
1 or 2 or 3 ...), in which case several thousand nodes 
may be treated. However, for irregular lattices and/or 
systems with variable coefficients, it may be necessary 
to store a separate set of coefficients for each node 
(d - m) , bringing the maximum number of nodes down to 
several hundred. 

(2) In order to obtain high accuracy with practicable lattices, 
it is frequently necessary to crowd nodes in regions where 



■12- 



a solution is known, either by analysis or previous numeri- 
cal experimentation, to have "large" third- or higher-order 
derivatives. In the neighborhood of a singularity, for exam- 
ple, the nodal density may need to be several orders of mag- 
nitude greater than in that of well-behaved regions. 

(3) In constructing lattices with variable nodal spacing, it 

§ 
is recommended that abrupt changes in spacing be avoided. 

This argues for the use of a functional mechanism for arriv- 
ing at nodal distributions. Such a mechanism may be speci- 
fied quite easily with the EPS command language, as is shown 
in Section 3.5. 

(4) In general, nodes should be arranged so that quadrilateral 
lattice elements are as nearly square as possible. Grid 
lines which intersect at angles less than 45° and lattice 
elements with aspect ratios greater than 5:1 should be avoid- 
ed in most applications, especially in sensitive regions. 

(5) In problems involving material interfaces or other conditions 
leading to differential equations with discontinuous coeffi- 
cients, it is generally more accurate and convenient to ar- 
range grid lines so that they coincide with lines of discon- 
tinuity. Similarly, nodes should be placed at all boundary 
points where boundary conditions change abruptly. 



+ 
See Parmelee [4], pp. 75ff. 

§ See Varga [5], p. 191. 
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(6) Other matters being equal, it is advisable to avoid lattice 
configurations which, because of the presence of diagonal 
boundary segments or superimposed nodes, contain triangular 
elements. Superimposed nodes should specifically be avoided 
at interior points of a domain, and should be used at bound- 
ary points only when boundary values, as opposed to gradi- 
ent or more complex conditions, are prescribed for all de- 
pendent variables. 

The above remarks are necessarily qualitative, for answers to specific 
questions, e.g., "How much accuracy can I expect from a lattice like 
'this'?", are available only for a few special classes of lattices, 
and then only for simple differential systems. In practical applica- 
tions it is generally necessary to determine solution accuracy empir- 
ically, for example, by treating for comparison a problem variation 
having a known solution or by modeling the actual problem with two or 
more lattices, each having a different number of nodes. In the latter 
approach, it may be possible to estimate errors roughly by noting the 
degree to which the numerical solution changes from lattice to lattice. 
In fact, if lattice modifications are handled in a systematic manner, 
the basic form of the lattice being held constant as more grid lines 
are introduced, it may be possible to obtain a greatly improved solu- 
tion by extrapolation. 



+ See Crandall [6], pp. 171-173 and pp. 269-271, and Parmelee 
[4], pp. 65ff. 
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3.3. Computer input: basic approach. As Section 3.1 may sug- 
gest, the process of selecting a lattice for EPS generally involves 
at least an implicit assignment of grid coordinates, it being conven- 
ient, in general, to introduce grid coordinates in order to test a lat- 
tice for admissibility. Grid coordinates are, in any event, essential 
to the procedure by which lattice information is communicated to the 
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Fig. 3. Lattice for a multiply -connected domain in (a) plane 
of independent variables x and, y, and, (b) i,,j -plane. 



computer, and should be assigned in a manner which is natural and con- 
venient to the problem at hand. In essence, the grid coordinates form 
a discrete, curvilinear reference system, making it possible to view 
any lattice as a mapping from the i ,,j -plane of a regular, orthogonal 
grid bounded by horizontal, vertical, and +45° lines (see Figure 3). 
The computer input required to define a lattice closely parallels this 



It is permissible to use both positive and negative integers 
provided magnitudes are kept less than 2 27 (134,217,728). 
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conceptualization, for it is necessary to provide the machine with a 
description of the lattice as it appears "in i,J-space" and with a map- 
ping algorithm by means of which i,j -values can be converted into cor- 
responding values for a problem's actual independent variables. The 
specific commands needed for these operations are discussed below. 

3.4. Specifying lattice structure in 7),j'-space. Lattice struc- 
ture in i, j -space is specified by defining the path of each bounding 
contour (only one in singly-connected domains) in terms of its i,j- 
coordinates. Specifically, it is necessary to type a command sequence 
of the form 

append i ± ,j 1 , i 2 ,0 2' •••> ^N'J'n to curvename 
close auTvenarne 

for each contour curvenarne , where grid coordinates i-pjj* ^2 "72' •••> 

i-n,jj\j represent the contour's vertices in the i,j'-plane, listed in 

t 
either clockwise or counterclockwise order. For example, for the lat- 
tice of Figure 3, the user could type 

append 0,0, 0,5, 2,7, 2,l6, 11, 16, 11,0 to alpha 

c 1 ose a 1 pha 

append 5,2, 5,4, 8,4, 8,2 to beta 

close beta $ 

The computer would respond by printing 



t 
Note that vertices in the ^,J -plane may or may not coincide with 

boundary discontinuities in the plane of the independent variables. 



■16- 



DEFINITION OF NEW CURVE 'ALPHA' HAS BEEN COMPLETED. 
'ALPHA' HAS BEEN CLOSED. 

DEFINITION OF NEW CURVE 'BETA' HAS BEEN COMPLETED. 
'BETA' HAS BEEN CLOSED. 

The use of "alpha" and "beta" as contour identifiers is, of course, 
completely arbitrary; in general, contours are named like algebraic 
parameters (see Section 2), it being essential only that each contour 
be given a unique identifier. 

3.5. Specifying a mapping algorithm. With EPS, the symbols "x" 
and "y" are always used to denote a problem's independent variables, 
even when these variables do not represent Cartesian space coordinates 
(see Section 4), and the procedure for specifying a mapping algorithm 
always involves at least a command of the form 

define x=f, y=g 

where f and g are expressions which either explicitly or implicitly 
depend on grid coordinates "i" and "j". In the simplest cases, i,,j- 
dependence can be shown explicitly, so that a single command of the 
form shown above is all that is necessary to establish the required 
mapping. For example, to define a regular, orthogonal lattice where 
the nodes are spaced 0.1 unit apart in the x-direction and 0.15 unit 
apart in the y -direction, and where the node with grid coordinates 
(0,0) corresponds to (x,y) = (0.3,-0.9), the user could type 

define x=i/10+.3, y=.15--(j-6) $ 
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In other cases it is convenient or necessary to show i , j -dependence in 
a more implicit manner, through the use of auxiliary parameters . For 
instance, a lattice representing the discrete counterpart of a cylindri- 
cal grid with spacing hv = 0.5 and A9 = 15° could be defined as follows 

define x=rad i us-cos (t heta) , y=rad i us»s i n (theta) , 

radius=.5- c i , theta=j- ; p i /12 
set pi= z t-atan (l) $ 

Note that the def i ne-command is used to indicate functional relation- 
s- 
ships while the set-command is used to specify constants. 

Both of the examples above represent very simple situations, thus 
are not completely germaine to many practical applications, where it is 
often necessary to introduce lattice irregularities in order to achieve 
accurate results. Irregularities may, of course, be of varying com- 
plexity, a particularly simple case being one in which the grid lines 
are straight and parallel to the x- and y-axes, but unevenly spaced 
(see Figure 2b). In instances where the grid-line spacing cannot be 
expressed algebraically in terms of "i" or "j", such configurations 
can be defined by associating with "x" and "y" vectors which give the 
desired spacing, as in the following sequence: 

define x=xx ( i ) , y=yy ( j ) 

set xx(0)=.53, xx(l)=.77, xx(2)=1.12, ..., 
yy(0)=.20, yy(l)=.42, yy(2)=.75, ■•• $ 



'Under normal circumstances the computer prints no explicit re- 
sponse in connection with define- and set-commands, but indicates com- 
pletion of these commands simply by printing another "PROCEED:". This 
is also true of several other commands and can be assumed to be the 
case in general when the text makes no mention of a response. 
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Similarly, it is possible to associate "x" and "y" with matrices so 
that the coordinates of each node can be specified individually: 

define x=xx(i,j), y=yy(i,j) 

set xx(0,0)-5.2, yy(0,0)-3.8, xx(l,0)=6.3, yy(l,0)=3.6, 
xx(2,0)-7.3, yy(2,0)=3.5, ... $ 

The typing involved in specifying the components of vectors and matrices 
can be reduced considerably in comparison to that shown above through 
use of the EPS macro facility (see especially Section 9.4). However, 
representing "x" and "y" with matrices is demanding at best. The user 
tempted to use such a procedure is urged therefore to note that it is 
frequently possible to have the computer generate internal node coordi- 
nates automatically by interpolating between values specified along the 
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Fig. 4. Lattice with linearly inter- 
polated nodal distribution. 



lattice boundary. For example, consider linear interpolation, as illus- 
trated by the lattice of Figure 4. This lattice could be specified by 
a sequence such as 
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define x=xO ( i ) +j * (x5 ( i ) -xO ( i ) ) /5 , 
y=yO(i)+j*(y5(i)-yO(i))/5 
set xO(0)=1.3, yO(0)=1.0, x0(l)=1.8, yO(l)=1.0, 
xO(2)=2.A, yO(2)=0.9, ... $ 

this being the discrete counterpart of a two-dimensional ruled sheet 
transformation, applicable to domains which have two or more non-adja- 
cent straight sides. A more general interpolation scheme, applicable 
to arbitrary domains, can easily be inferred from the surface equations 
discussed by Coons [7] . 

As a final point it should be noted that when a lattice assumes 
distinctly different forms in different regions, it is possible and 
sometimes convenient to introduce branches in the mapping algorithm so 
that mechanisms corresponding to the different regions can be specified 
independently. In general this is accomplished by using relational and 
Boolean operators (see Table 3, Section 16) to form expressions which 
take on different discrete values for different ranges of "i" and "j". 
Suppose, for example, that an orthogonal lattice is desired in which 



x = .5 i, y 



.25 j for j < 8, 

.5 (j-8) +2 for 8 < j < 20, 

.75 (j-20) + 8 for j > 20. 



This lattice can be defined in terms of the relational operator "grt", 
which in an expression of the form "a grt 3" produces the value 1 when 
a > 3 and the value otherwise. Specifically, we can type 

define x=.5*i, y=yy((j grt 8)+(j grt 20)), yy(0)= 
•25*j, yy(l)=.5*(j-8)+2, yy (2)=.75*(j-20)+8 $ 
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I'his example, chosen for convenience, is not meant to encourage the use 
of lattice.- uhic. ch.ii^e character abruptly (see recommendations of 
Section .1.2). !he technique it demonstrates is perhaps most useful in 
situation? whore: ;;?""■ :: r :: lattice can be defined in a simple, alge- 
braic mri.-ir.' 1 :" , while .c"'u.-c earts require tabular definitions. This is 
the case, for example, in tee lattice of Figure 2a. 



4. Equation specification 

4.1. Admissible forms. The difference equation generator em- 
ployed by EPS models linear systems of partial differential equations 
given in terms of dependent variables u , u~, ..., u and independent 
variables x and y as follows: 



l { fe (« fef ^l + b kl ^1 + % Z u t ) 



1 = 1 " dx dy 

Z + e 

dx ""*" ciy 



(1) 



;— (d, Z + e,p Z+f LD Up)+a L pUp] = h, 
dy kZ t: — kZ -z — - kZ Z "kl Z> k 



[k = 1,2,. ..,n] 



Here the a, », b,p, ..., g,n and h, represent constants or functions of 
x and y which, along with the number n, must be specified by the user 
to indicate the particular system of interest. Note, for example, that 
Eq. (1) reduces to Laplace's equation in Cartesian coordinates, 



32 f + S - 0, (2, 



when 



dx z dy' 



u 1 -> 9, 
n - 1, 

a ll = e ll = lj 



(3) 



b = a = d = f = a = h = 0. 
11 11 11 J ll «11 1 



It should be emphasized that the independent variables x and y in Eq. 
(1) may denote position in any right-handed coordinate system. Thus 
Eq. (1) corresponds to Laplace's equation in cylindrical coordinates, 
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d r 8$> . 1 ^ 2 



+ ^2 = 0, (4) 



when 

x -> p, y -> 8 , u, -> $, 
n = 1, 

*ll = c ll = d ll = ^11 = ^11 = ^1 = °' 



(5) 



Such a particularization is admissible, of course, since the a^, fc^, 
etc. may vary with position; as a result, the user is free to employ 
whatever coordinate system best suits his problem. More importantly, 
he may deal with physical media whose properties change from point to 
point . 

It is important to note that while the EPS difference equation 
generator models linear differential systems, solutions to nonlinear 
systems may frequently be obtained without difficulty. In essence this 
is done by having the computer produce and solve a sequence of linear 
problems whose spatially varying coefficients represent successively 
improving approximations to the solution-dependent coefficients for 
the nonlinear problem of interest. In this section we are, in fact, 
concerned primarily with the solution of linear equations, it being 
desired at this point only to indicate the feasibility and general char- 
acter of the procedures needed for nonlinear applications. Further de- 
tails and examples of such applications are given in Section 10. 

4.2. Specification of parameters. The quantities n, a^, b^, 
etc. which indicate the particularization of Eq. (1) of interest are, 
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like the lattice parameters "x" and "y", specified through use of the 
set- and def i ne- commands . As mentioned earlier, "set" is used for the 
specification of constants while "define" is used to specify functional 
relationships, e.g., to show how equation coefficients change with posi- 
tion. Before specific examples are discussed, the following remarks 
should be noted. 

(1) When n > 1, equations should be ordered so that in the k 

equation u, is the "dominating" dependent variable. Specifi- 



cally, it is mandatory that u, appear in the k equation, 



k 

k 

i.e., that at least one of the a,,, b,,, ..., g<< be non-zero, 
and it is desirable that the following inequality be satisfied: 

Kfe + e fefcl - Wl + e kl\ lk,l= l,2,...,n] 

(2) As Section 3.5 implies, EPS uses standard algebraic compiler 
notation for subscripted variables. Thus a„_ is typed "a (2,3)" 
and h 2 as "h(2)". 

(3) When all components of a subscripted variable are identical, 
subscript notation may be dropped and all components specified 
at once. For example, if coefficient a, „ = for all k and I, 
one may simply "set c = 0". This means, in particular, that 
when n = 1, subscripts may be omitted for all equation parame- 
ters . 

(4) Spatially varying equation parameters can be specified in terms 
of grid coordinates "i" and "j" and/or in terms of any quanti- 
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ties which depend on "i" and "j", e.g., "x" and "y". More- 
over, the functional behavior of a parameter can be given 
in algebraic, tabular, or mixed algebraic and tabular form. 
In any case, it is important to recognize that a discretely 
varying quantity such as "x" will be modeled as a continuous 
variable only when it operated on by the built-in interpola- 
tion function "fit", as in "fit(x)". We consider two common 
situations : 

(a) Algebraic expressions in x and y. When one or more 
equation parameters are known algebraically in terms 
of the independent variables x and y, auxiliary vari- 
ables representing "fit(x)" and "fit(y)" should be de- 
fined and used in place of x and y when the parameters 
in question are specified. For example, if 622 = ^'^ 



+ x' and h. = x + Jy, type 



define fx=fit(x), fy=fit(y), e(2,2)= 
3.5+fx power .8, h(3)=fx+sqrt(fy) 

(b) Tabular representations. When an equation parameter 
a = a(z), where 3 is a function available only in dis- 
crete form, say as a matrix "z(i,j) M with components 
which give s's values at the lattice nodes, then an 
auxiliary variable representing "f i t (z(i , j) )" should 
be defined and used in the specification of a: 

define fz=f i t (z(i , j)) , a=aTfz; 

Here "z(i,j) M could represent information determined 
by experimental measurement and specified component- 
by-component to the computer, or it could represent 
discrete information calculated by the computer, as 
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is the case in the solution of nonlinear equations 
(see Section 10.2). 

(5) Discontinuous equation information is specified by using the 

relational and Boolean operators in a manner similar to that 

introduced in the discussion of lattice definition procedures 

in Section 3.5. Suppose, for example, that h = 5 for x < 

7.2 and h = for x ^ 7.2. This behavior could be specified 

in terms of the relational operator "les" (read "less than") 

by typing 

define fx=fit(x), h(2)=hh(fx les 7.2) 
set hh(0)=0, hh(l)=5 

or, alternatively, by typing 

define fx=fit(x), h(2)=5*(fx les 1.2) 

As it happens, specifying discontinuous parameters in terms 
of "x" and "y" results in a certain amount of smoothing when 
lines of discontinuity do not coincide with grid lines. More- 
over, for arbitrarily curved lines of discontinuity it is dif- 
ficult to express discontinuous behavior in the manner shown 
above. These disadvantages may be avoided if, as recommended 
earlier, grid lines are made to coincide with lines of dis- 
continuity. When this is the case, it is possible to express 
discontinuous parameters in terms of "i", "j", and a built-in 
(argumentless) function "octant", which, when difference equa- 
tions are being formed, indicates lattice sector identity. 
Specifically, "octant" takes on integral values between and 
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Fig. 5. Numbering of lattice sectors. 

7, as shown in Figure 5, and may be used to differentiate 
between one side of a grid line and the other. For example, 
if h- = 5 for lattice elements below j = 20 (in -£,j-space), 
and h = for elements above j = 20, one can type 

define h(2)-hh((J les 20)+(j les 2 )), 

hh(l)=5* (octant les 4) 
set hh(0)=0, hh(2)=5 

This notation is tedious, but quite flexible and easily mas- 
tered with a little practice. 

4.3. Examples. Shown below are input sequences corresponding 
to four special cases of governing system (1). Additional examples 
are given in Section 10.2. 



mand 



A. Laplace's equation in Cartesian coordinates. The single corn- 



set a=e=n=l, b=c=d=f=g=h=0 $ 
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can be used to specify Laplace's equation in Cartesian coordinates. 
(See Eqs. (2] and (3).) 

B. Laplace ' s equation in cylindrical coordinates . The command 

sequence 

set n=l, b=c=d=f=g=h=0 

define fr=fit(x), a=fr, e=l/fr $ 

can be used to specify Laplace's equation in cylindrical coordinates. 
(N.B., the independent variables are r and 9, not r and z; see Eqs. 
(4) and (5).) 

C. Plate equation in Cartesian coordinates. The deflection 

w = w(x,y) of a flat plate with variable stiffness s = s(x,y) subjected 
to a pressure field p = p(x,y), where x and y are Cartesian coordinates, 
satisfies 

div (s grad V 2 w) = p 

for sufficiently small p/s. This fourth-order system may be rewritten 
in the form of Eq. (1) on introduction of an auxiliary variable = V w, 
as follows: 

div (fl grad $) = p 

$ - V 2 w = 



With $ and w associated with u and u 9 , respectively, the command se- 
quence to prescribe this system becomes 



set n=2, a(l,2)=a(2,l)=b=c=d=e(l,2)=e(2,l)=f=g(l,l)= 
g(l,2)=g(2,2)=h(2)=0, a (2 ,2)=e (2 ,2)=-l, g(2,l)=l 
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define fx=fit(x), fy=fit(y), a(l,l>=e(l,l) a! sff^fy) , 
h(l)= p(fx,fy) $ 

(Note that s(fx. 3 fy) and p(fx y fy) represent expressions to be filled in 
by the user.) 

D. Plane stress equations expressed in terms of displacement 
variables 3 Cartesian coordinates. The Cartesian displacements u = 
u(x,y) and v = v(x,y) for a plane elastic continuum with constant mate- 
rial properties E and v satisfy 

tel 1=^ 2 ( to + V 9^ ^ + ^ 27TkO ( ty to )J 

(6) 

to*- 2(l^y ( 9^ + to ^ + 9^ l^J 2 ( V to + 3y } J " ° 

when no stresses normal to the x,z/-plane and no body forces are present. 
With u and y associated with «, and «2» respectively, the command to 
specify this system becomes 

set n=2, c=f»g=h=a(l,2)=a(2,l)»b(l,l)='b(2,2)«d(l > l) = 
d(2,2)=e(l,2)=e(2,l)=0, a(l,l)«e(2,2> £/(l-v*v), 
b(l,2)=d(2,l)= v*a(l,l), a(2,2)=b(2,l)=e(l,l)= 
d(l,2)- E/(2*(l+v)) $ 

XHere E and v represent numerical values to be filled in by the user.) 
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5. Specification of boundary conditions 

5.1. Admissible forms. The EPS difference equation generator 
admits boundary conditions for governing system (1) which designate 
values for, or linear relationships between, the dependent variables 
U/> and certain first-order functions F, [k = 1,2,..., n] of the w». 
The F, may not be elected by the user, but are restricted to a prede- 
termined form which is natural and convenient for most applications 
(see Section 5.5). Specifically, the F, are given in terms of the co- 
efficients a^n, t>uo, ■ • • , fuo of governing system (1) and the local 
boundary slopes dx/a& and ay/ai>, arc length 4 being positive when the 
boundary is traced so that the domain lies to the right, as follows: 

' t ■ I ( - *** £t ♦ »u £* + "u «*> £ 

£=1 ax ay ,_.. 

+ (d kt ^+e kl ^l + f bl u t ) |f ) 

ax ay 

[k - 1,2 n) 

Note that for Laplace's equation in Cartesian coordinates (see Eqs. 
(2) and (3)), definition (7) reduces to 



au. ,. au, ~ du- «. 

_ _ 1 3i£ 1 3x m 1 ^ 3£ 

1 = " ax 34 ay a6 " an an 



where n represents the direction of the local outward normal. Thus, 
for this special case , boundary conditions must assume the form 



a(x,y) || + Hx,y) $ - y(x,y) (8) 
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where either a(x ,y) or $(x,y) may become zero. More generally, for 
each boundary segment m, boundary conditions for governing system (1) 
must take the form 



n 
J/ ?\dm F l + %Zm u l> = r km (9) 



[fe = 1,2, 



where the p,p , q,p , and v, are constants or functions of position 
to be specified by the user. The specification of boundary conditions 
therefore entails two distinct operations: 



(1) assignment of unique segment numbers to portions of the 
boundary having different types of boundary conditions; 

(2) specification of parameters p,n , q<p , and r, for each 
segment number m to indicate the desired condition. 

The specific input conventions by which these operations are effected 
are discussed below. 

5.2. Assigning segment numbers. The boundary of the domain 
should be divided into segments, these corresponding to arcs over which 
the boundary conditions assume constant form, i.e., continuous behavior 
in terms of the p's and q ' s of Eq. (9). Each segment should then be 
assigned a number, say m, and identified to the computer in terms of 



t 
Other forms, including nonlinear relationships between the F^ 

and up, can be accomodated by iteration. However, as explained in sec 

5.5, other forms seldom arise in physical applications, even when the 

governing equations are nonlinear. 
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the grid coordinates i-,,,3-, and i j of its end points plus the coordi- 
nates ^ q of some intermediate point, as follows: 

impose m along i } ,j } , ? : 2 " 7 '2' i 3" ] '3 
The following remarks should be noted: 

(1) Segment numbers may be arbitrary integers, but, to conserve 
memory, should be assigned in sequential fashion. Generally 
it is convenient to use the numbers 0,1,2,... (The number 

is automatically assigned to segments not explicitly men- 
tioned in impose-commands . } 

(2) A given segment number may be used more than once if the same 
boundary condition applies to more than one portion of the 
boundary. 

(3) The i,,j, may refer to arbitrary boundary nodes; they need 
not refer to the contour vertices discussed in Section 3.4. 

(4) Intermediate coordinates i ,j' ? are omitted in the case of 
a segment which joins neighboring boundary nodes. 

These points are illustrated by the following command sequence, which 
could be used to assign the segment numbers shown in Figure 6. 

impose 1 along 0,0, 1,0, 9,0 

impose 2 along 0,0, 0,1, 6,7 

impose 3 along 6,7, 7,7 

impose 2 along 9,0, 9,3, 7,7 

impose k along 3,3, 3,5, 3,3 
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Fig. 6. Sample lattice and segment numbers (in brackets) . 
Ends of segments are indicated by large dots • . 



5.3. Specification of parameters. The parameters Fy^t Qhpm* 
and r, of Eq. (9) are specified in exactly the same manner as the 
parameters a,p, J, », ..., h, of governing system (1) and, in particu- 
lar, are subject to the same input conventions as described in points 

x 
(2) -(5) of Section 4.2. Thus, for example, to specify q l2 % = 1 + e , 

type 

define fx«fit(x), q (1,2,3 )=l+exp(fx) 

and so forth. Users with straightforward applications, especially 
problems with only one dependent variable (n = 1), may feel free at 
this point to skip to Section 5.4. For more advanced applications, 
the more detailed remarks below may be relevant. 
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(1) Because of the manner in which EPS models boundary condi- 
tions, greatest accuracy is achieved when the coefficients 
Phlm are cons,tant; along each- boundary segment w». Thus, if 
possible, any space dependence which might naturally arise 
in the Vugm should be eliminated by division. For example, 
with Laplace's equation in Cartesian coordinates, where the 
boundary condition is as shown in Eq. (8) with a(x,y) f 0, 
use 

Clearly such a procedure is always possible fo« n = 1, and 
it is also possible for n > 1 so long as boundary conditions 
do not involve linear combinations of the £*> 

(2) When n > 1, equations that express boundary conditions should 
be ordered, if possible, so that for each boundary segment m 
and each equation k at least one of the coefficients Puuyy. and 
^kkm * s non-zera< That is, if possible, the k equation 
should involve either F. or u^. or both. If for any m such 
an ordering is impossible, make certain that the boundary 
conditions are well posed, i.e., that they are self -consist- 
ent and sufficient to determine all of the -Uu uniquely, then 
proceed with an ordering which maximizes the number of equa- 
tions with non-zero qkhm' 

(3) At a node 8 which separates boundary segments with different 
segment numbers m.. and m~, the coefficients puffo for each 
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equation k must satisfy at least one of the following con- 
ditions : 

(P klmi h = ° ^ = 1,2,..., n] 

or 



2 



or 



Thus, for example, in boundary-value problems for Laplace's 
equation, boundary conditions may switch abruptly from Neu- 
mann conditions (v^^ 4 0, a-,-, =0) to Dirichlet condi- 
^llm ' '11m 

tions (p.... = 0, a-,-, 4 0) , but not from "Neumann condi- 
VK llm ' Y llm ^' 

tions to Neumann conditions" unless p.... = p,, at the 

L 11m i ll>n 

point of discontinuity. Again, for n = 1 there need be no 
conflict; in fact, the restriction cited here is satisfied 
automatically if boundary conditions are normalized so that 

the p,, assume only values of unity or zero. However, 

-llm J ; 

when n > 1, conflicts may arise in rare cases. (If not an- 
ticipated by the user, they will be detected and reported 
by the computer at the time an attempt is made to generate 
difference equations.) In such cases, it is necessary to 
force the computer to employ the same segment number, say 
m , on both sides of the node in conflict. This is done 
by using a special form of the impose-command 

impose m at £ R ,Jt> 
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which effectively extends segment m 1 half-way to 8's neigh- 
bor on the "m^-side." In most instances, this effect may 
be counteracted by proper arrangement of nodes, but, at 
worst, it merely introduces an additional discretization 
error . 

(4) Boundary condition parameters may be expressed in terms of 
boundary slopes 9o;/34 and 3z//34 through use of the built-in 
(argument less) functions "xs" and "y s "» respectively. Thus 
if r 12 - (dx/d&)/y, type 

define fy=fit(y), r(l,2)=xs/fy 

5.4. Example. Consider the solution of Laplace's equation in 
terms of the lattice and boundary conditions illustrated in Figure 7. 
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Fig. 7. Sample lattice and boundary con- 
ditions for Laplace 's equation. 
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The boundary conditions shown here can be specified as follows: 

impose 1 along 0,0, 0,2, 0,4 
impose 2 along 3,4, 4,4, 6,4 
impose 3 along 0,4, 2,4, 3,4 
impose 3 along 0,0, 5,0, 6,4 

set p(l,l,l)=p(l,l,2)=q(l,l,3)=r(l,l)=r(l,3)=0, 
q(l,l,l)-q(l,l,2)-p(l,l,3)-l, r(l,2)=l00 $ 

Significant reductions in the amount of typing required to specify 
boundary conditions can be realized by means of appropriate macros. 
See Sections 9.1 and 9.7. 

5.5. Supplementary remarks. For users with applications involv- 
ing governing expressions other than Laplace's equation, an explana- 
tion regarding the fixed form of the functions F, shown in Eq. (7) 
may be desirable. For this purpose, it is helpful to consider the 
special case in which the independent variables x and y represent Car- 
tesian coordinates. This makes it possible to express the local unit 
outward normal to the boundary in terms of boundary slopes 9:r/94 and 
9y/9<6 and unit vectors i and j, corresponding to the x- and ^-direc- 
tion, respectively, as follows 

n - ' x af + J 9T 

and suggests that the F. be written as vector dot products 

\ - F fe ' n (10) 

where the vectors F. are defined by 
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£=1 ax ay 

n ' (11) 

t-1 ox dy 

[k =» 1,2, ... ,w] 



Definition (11), in turn, suggests a more succinct statement of govern- 
ing system (1), viz., 

_ n 
div F fe + I gki u i " h k [fe -■ 1,2,...,*] (1*) 

This will be seen to resemble closely expressions commonly used in the 
representation of physical laws of continuity and equilibrium. Thus, 
in physical applications the vectors Fl are associated with mass flow 
vectors, current density vectors, stress dyadics, and so forth, while 
the related functions F, represent boundary fluxes, currents, mechani- 
cal tractions, and the like. From this we may draw two important con- 
clusions: 

(1) the F, are, at least in physical problems with Cartesian 
coordinates, exactly the quantities with which we normally 
prefer to deal, and 

(2) the physical significance of the F± can frequently be de- 
duced without close scruntiny of the formal definition 
given in Eq . (7) . 

In fact, it will be found that these conclusions are valid, with minor 
qualifications, for any orthogonal curvilinear coordinate system. 
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Moreover, they are generally valid regardless of whether the govern- 
ing system being treated is linear or nonlinear. Nonlinearities in 
physical problems normally arise in the constitutive laws which relate 
vectors F, (or their counterparts in non-Cartesian reference frames) 
to the dependent variables u„, not in the equations which govern the 
actual behavior of the F, . 

A final example, which illustrates point (2) above, can be given 
in terms of the plane stress equations introduced in Section 4.3 (Exam- 
ple D, Eq. (6)). These equations merely state the equilibrium condi- 
tions 

9a 3t 

dx dy 

(12) 

9t do 

__xy + _yy = Q 

dx dy 

for two-dimensional Hookean elasticity. Comparing Eq. (12) to Eqs . 
(10) and (1'), we can deduce immediately 

, - - . - boundary traction 

1 v xx J yx in ^-direction 

. - - . - boundary traction 
F = ( i t -Hie ) * yi — 

2 xy J yy in y -direction 

Accordingly, we see that the general form for boundary conditions for 
problems with two dependent variables, i.e., 

Pllm F l + Pl2m F 2 + q -llm U l + ^12m U 2 = V m 
Pllm F l + P22m F 2 + q 21m U l + < *22m U 2 = r 2m 
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becomes, in the case of plane elasticity problems, a pair of linear 
equations relating the boundary tractions and displacements. These 
equations may be specialized so as to designate simple traction and 
displacement conditions, or may be used in a more general fashion to 
treat boundaries on elastic foundations, and so forth. 
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6. Getting a solution 

6.1. Generation of finite-difference equations. Once a lattice 
and differential system have been specified as indicated in Sections 3, 
4, and 5, or revised as indicated in Sections 8 and 10, finite-differ- 
ence equations can be obtained simply by typing 

form $ 

Assuming there are no detectable inconsistencies in the user's data and 
difference equations are being obtained for the first time, this command 
will result in an immediate response of the form 

POINT TALLY IS N. 

SPACE FOR SOLUTION MATRIX HAS BEEN ALLOTTED AND ZEROED. 

where N is the total number of nodes in the lattice. When equations are 
being regenerated after a problem revision, the response may be abbrevi- 
ated or omitted entirely. In any event, a pause will ensue which may, 
in fact, be imperceptible or last for several minutes, depending on 
problem complexity and computer load. Eventually a PROCEED-message will 
be printed, at which point it may be assumed that all difference equa- 
tions have been formed. 

6.2. Solution of finite-difference equations. EPS solves dif- 
ference equations generated by means of the form-command by successive 

+ 
relaxation, using a relaxation factor "omega", a tolerance parameter 



4- 

See, for example, Young [8], sec. 11.5. 
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"delta", and a maximum iteration index "limit" specified by the user. 
Relaxation factor "omega" has the usual significance, thus may be as- 
signed a value between and 1 for underrelaxation or a value greater 
than 1 for overrelaxation , according to the situation at hand (see 
below). Relaxation, once initiated, continues until either 

(1) the absolute changes introduced in all components of the 
solution matrix during a given sweep of the lattice come 
to be less than or equal to "delta", or 

(2) the number of sweeps reaches "limit". 

After this, it is possible to request further relaxation, perhaps with 
different values for "omega", "delta", and "limit", or to proceed imme- 
diately with requests for output. It is also possible to intersperse 
requests for output with requests for further relaxation in order, for 
example, to monitor the convergence of individual components of the 
solution 

Parameters "omega", "delta", and "limit" are specified by means 
of the set-command, as in 

set omega=1.6, delta=2e-5, limit=200 

while relaxation is initiated or resumed by typing 
relax $ 

The rel ax-command always results, after one or the other of the termi- 
nation conditions has been satisfied, in a response of the form 

RELAXATION TERMINATED AFTER N PASSES. MAX SOLN CHANGE: C. 
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where N is the number of sweeps of the lattice actually performed and 
C is the ™* X . |Au(k,£,j)| for the Nth sweep. 1 " Note that either N will 
equal "limit" or C will be less than or equal to "delta". 

Selection of "omega". In general, there is for any set of dif- 
ference equations a maximum relaxation factor u> max beyond which succes- 
sive relaxation diverges. In many common applications U) max is 2, while 
convergence is most rapid for some optimum relaxation factor between 1 
and 2. Thus 1.5 is often a reasonable first choice for "omega". For 
lattices with many nodes (say more than 300) or in cases where a number 
of related problems are being treated, it may be worthwhile to attempt 
to estimate the optimum relaxation factor by experiment. Normally this 
can be done by setting "limit" to some relatively small number and ob- 
serving the rate of decrease of the "MAX SOLN CHANGE" for different 
values of "omega". § Such a procedure may also be required in certain 
cases where divergence occurs for values of "omega" in the normal range, 
It has been found, for example, that u) is less than 2 for certain 
plate problems with clamped-edge boundary conditions (see Section 4.3, 
Example C) and considerably less than 1 for certain nonlinear systems 
arising in hydrodynamics. 



Here, as in actual computer input, components of the solution 
matrix are denoted by "u {k,i,3)" t meaning the current approximation 
for dependent variable u, at node (i,{j). 

Convergence rates are often dramatically dependent on the relax- 
ation factor. However, for lattices with no more than a few hundred 
nodes, the actual computation time required for relaxation with any 
"reasonable" factor is typically much less than that required for other 
operations involved in obtaining a solution, e.g., the generation of 
difference equations. In such cases there is little to be gained on 
a percentage basis by seeking the optimum factor. 
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Selection of "delta". It is emphasized that tolerance parameter 
"delta" can be used only to control errors in the discrete solution 
u(k,i,j) that are associated with incomplete relaxation; its value has 
no bearing whatever on the discretization errors incurred by replacing 
a continuous domain with a discrete model. Moreover, "delta" is only 
indirectly related to relaxation error, its strict function being to 
specify a maximum allowable value for the "MAX SOLN CHANGE" computed 
during a given sweep of the lattice. Normally the maximum change oc- 
curs in a component of the solution matrix which itself has a relative- 
ly large absolute value. Thus assigning "delta" some value d means 
loosely that the relaxation should be terminated if and when the larg- 
est components of the solution matrix are changing by amounts no great- 
er than d. It does not mean, however, that when the relaxation is 
terminated on this basis, the largest components will be "good to ±d" , 
for convergence to the true discrete solution may occur in steps which 
are much smaller than the actual error. For this reason, "delta" should 
generally be given a value which is one or two orders of magnitude less 
than the relaxation error considered tolerable in the largest solution 
components. Note that convergence is normally fairly uniform with 



+ 
This is a rule of thumb which should be appropriate for most 

applications. However, more conservative, i.e., smaller, values for 
"delta" may be advisable if convergence is abnormally slow, as when 
difference equations are ill-conditioned and/or it is necessary to 
resort to small relaxation factors, e.g., less than 0.5. The ratio 
of relaxation step size to actual error is, of course, the critical 
variable in all cases. This can be estimated, if doubts arise, by 
observing the rate of decrease of the "MAX SOLN CHANGE" as the relax- 
ation progresses. For example, if 80 sweeps are required to reduce 
this number by an order of magnitude, it is reasonable to assume that 
the reported solution changes are about 1/8 the size of the actual 
errors due to incomplete relaxation. 
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respect to all components of the solution, i.e., fractional rates of 
change in the smaller solution components can be expected to be about 
the same as those in the larger components. 

Selection of "limit". Parameter "limit" has two intended func- 
tions: (1) to prevent a nonconvergent relaxation from continuing in- 
definitely, and (2] to permit close monitoring of the relaxation, as 
is necessary in estimating optimum values for "omega", etc. Here it 
suffices to say that fully relaxed solutions, i.e., solutions for which 
further improvement is impossible because of the effects of computa- 
tional round-off error, can be obtained in typical applications with 
a number of sweeps which is about the same as the number of nodes in 
the lattice. As a rule, then, "limit" should not be given a value 
much larger than the number of nodes, although this may mean in some 
cases that additional re 1 ax- commands are needed to reach full conver- 
gence . 



If full convergence is sought, it will be observed that, after 
decreasing steadily for a time, fractional solution changes (and hence 
reported values for the "MAX SOLN CHANGE") ultimately reach a lower 
limit beyond which further relaxation is fruitless. In some cases 
this limit is about the same as the round-off error size intrinsic to 
the computer, i.e., about 10 -8 on the IBM 7094. However, it is not 
uncommon for the limit to be larger, e.g., 1CT 5 or 10 _b , especially 
with lattices of several hundred nodes. 
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7. Output 

EPS provides both printed and graphical output representing the 
following quantities: 

(1) solution components u(k ,i ,j) ; 

(2) numerically evaluated first and second derivatives ux(k,i,j), 
uy(fc,i,j), uxx(k,i,j), uxy(k,ij), uyx{k,i,j), and uyy (k ,i ,j ) ; 

(3) numerically evaluated boundary fluxes f]ux(k,i,j) correspond- 
ing to the F[ z of Eq. (7); and 

(4) arbitrary functions of any of the above as well as of other 
quantities, not specifically related to a finite-difference 
solution, which the user has defined. 

The particular commands by which output is obtained are described in the 
two sections which follow. 

7.1. Printed output. Small quantities of printed output are ob- 
tained by using commands of the form 

print e 1 , e 2 > • • • » e N 

where the eu are algebraic expressions whose values are to be printed in 

a standard scientific (!•-) format. For example, to determine current 

numerical approximations for w 2 ant ^ (S^/Sy) / i^u^/dx) at node (7,9), 
type 

print u(2,7,9), uy (3,7,9)/ux(3,7,9) $ 
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The computer will respond with a line of the form 

±X.XXXXXXX£±XX ±X.XXXXXXXE±XX 

assuming that the indicated quantities are defined (see Sections 14 and 

15). 

Larger quantities of output, e.g., tables representing solution 
properties over a wide range of i,j -values, can be obtained by embedding 
a print-command in one or more "do-loops". In general, an input sequence 

4* 

of the form 

do(s for v=v^ step 6v until v 2 ) 

can be used to cause repetition of any command sequence s while some 
parameter v is incremented by Sv from v ± to v 2 (inclusively) . Thus to 
cause printing of u 2 , Zu 2 /2x and du 2 /Zy along the grid line j = 5, for 
i = 1,2, ...,20, type 

do(print u(2,i,5), ux(2,i,5), uy(2,i,5) for 1=1 
step 1 until 20) $ 

This will result in the printing of 20 lines of output, each containing 
three numbers in the same format as shown in the preceeding paragraph. 
Similarly, to cause printing of u 1 and u 2 along the grid lines j = 2,4, 
6, 8,..., 20, for i - 8, 9,..., 11, type 



+ This input form does not represent a regular EPS command, but is 
transformed automatically into an appropriate sequence of regular com- 
mands in accordance with a macro definition which has been prespecified 
in both versions of EPS. See sec. 9.3. 
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do(do(print u(l,i,j), u(2,!,j) for i=8 step 1 until 11) 
for j=2 step 2 until 20) $ 

In the last two examples, no output of indices *'i" and "j" is requested, 
thus none is produced. For small tables of output this may be appropri- 
ate, but in other cases it is convenient to include "looping indices" in 
lists of items to be printed. As a completely independent matter, it is 
recommended that when complicated functions are being printed, they be 

defined beforehand as auxiliary parameters, so the auxiliary parameter 

t ... 

names can be used in actual print -requests. For example, if in solving 

a plane stress problem (see Eqs. (6) and (12)) numerical values for the 

stress components 0^, o y y, and t^ are desired along the grid line i = 

3, for 3 - 30, 29, 28,..., 5, type 

define sx= ff*(ux(l, i ,j)+v*uy(2,i ,j))/(l-v*v) , 
sy- F*(uy(2,r,j)+v*uxCl,ij))/(l-v*v), 
tau- ff*(ux(2,i,j)+uy(l,i,j))/(2*(l+v>) 

set i=3 

do(print j,sx,sy,tau for j-30 step -1 until 5)$ 

Parameters "sx", "sy", and "tau" could be used in any further requests 
without being redefined. 



The generation of large tables of printed output can consume 
considerable computation time, even as much as 0.5 second per line when 
complicated functions are involved. Using auxiliary parameter names in 
print-requests, rather than whole expressions, is one way of reducing 
output time requirements. Further savings can be realized by specifying 
and using looping macros which are less general, therefore more effi- 
cient, than the do-macro already provided. However, getting printed 
output with EFS is more expensive at best than, say, the user accustomed 
to batch processing with FORTRAN might expect. Discretion and selectiv- 
ity are therefore advisable. 
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While the standard format used by EPS for printed output is capa- 
ble of representing with full accuracy any number that EPS can produce, 
it is by no means an ideal format in all situations. For this reason, 
a special form of the print-command is recognized which allows the 
user to override the standard format with one of his own design. Spe- 
cifically, to indicate that output is to be printed according to a 
MAD- language format specifier /, type 

print e-p e 2 > •••> e /y using bed XfX 

f § 
where X is any single character not present in f. An example is 

provided by the sequence 

set i=j=7 

print x, y, u(2,7,7), flux(2,7,7) using bed Af9.4*/ $ 

which would cause to be printed a line of the form 

x.xxxx x.xxxx x.xxxx x.xxxx 

assuming all results are positive and less than 10. The standard output 



t 
MAD format conventions, described in [9], sec. 2.15, are virtu- 
ally identical to those of FORTRAN. N.B., with EPS, output lines may 
contain no more than 83 characters, format variables are not allowed, 
and carriage control characters are meaningless (the first characters 
of output lines are printed literally). 
§ 
Expressions of the form "bed XsX" , called literals , must not be 

split up over more than one physical input line. The computer inter- 
prets the character X, which, by definition, is the second character 
following reserved symbol "bed", as a logical quotation mark, and treats 
string s as a single (logical) word. Literals provide an escape mechan- 
ism which can be used not only in conjunction with format specifiers, 
but any time it is desired to have treated as a single word a character 
string which would normally be illegal or, at any rate, not interpreted 
as a single word. See sec. 15.1. 
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format allows for a maximum of five numbers per line, overflow being 
carried to additional lines when print- commands specifying more than 
five expressions are employed. Obviously an overriding format, speci- 
fying smaller field widths and more numbers per line, can be intro- 
duced if the user desires. 

7.2. Graphical output. When used in cpnjunction with the ESL 
display equipment, EPS (as represented by the disk file "eps-p saved") 
will accept a command of the form 

plot z 

where z is an expression giving the elevation at node 0,j) of a sur- 
face to be represented in graphical form. For example, to obtain a 
display of the finite-difference solution for dependent variable u. = 
u^x.y) , type 

plot u(l,i ,j)$ 

2 2 5 
Similarly, to see [(3u„/3ar) +(9u„/3z/) ]' , type 

plot sqrt(ux(2,i ,j)power 2+uy(2, i , j)power 2)$ 

Surfaces are represented as orthographic projections of the three- 
dimensional nets produced by moving the nodes of a finite-difference 
lattice out of the x,z/-plane by amounts proportional to argument z 
(see Figure 8). Projection planes can be varied continuously, i.e., 
a surface can be made to appear to rotate in space, by means of the 
plexiglass globe adjacent to the display unit. Rotating the globe 
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Fig. 8. Typical displays (photographically reversed) 
produced by means of the plot-command 



about one or more of its three axes causes the surface to rotate in 

t 
a corresponding manner. Releasing the globe stops the rotation. 

The user will find that the first display he requests appears 
with a "top-view" orientation, i.e., it shows the surface as seen by 
looking down the z-axis. Bach subsequent display, on the other hand, 
appears with the same orientation as that in which the previous dis- 
play was left unless an explicit command, namely, 

plot ni 1 S 

is given to cause all rotation angles to be reset to their initial 
values. The above command represents one of several special forms 
available for altering the normal operation of the display mechanism. 
Other forms can be used, for example, to suppress certain picture ele- 
ments, inhibit automatic scaling, etc. For details, see Section 15. 



In addition to the globe, the user will find near the display 
unit a control panel containing three large knobs. These knobs can 
be used to center the display or to change its magnification, although 
it is normally unnecessary to make such adjustments. 
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8. Data modification 

Previous sections of this document introduce five commands used 
to specify problem input. These commands can be grouped conveniently 
as follows: 

(1) the contour description commands, "append" and "close", 
used to define lattice structure in -£,j-space, 

(2) the parameter description commands, "set" and "define", 
used to specify various lattice, equation, boundary 
condition and relaxation parameters, and 

(3) the segment identification command, "impose", used to 
number boundary arcs having different types of bound- 
ary conditions. 

The present section describes how data associated with each of these 
command groups can be modified, either for purposes of error correc- 
tion or to introduce problem variations. 

Two general comments should be noted at the outset. First, it 
should be recognized that with EPS data can be modified in a very speci- 
fic manner. The system always maintains intact information which is not 
explicitly involved in a modification. Thus it is seldom necessary to 
"start over from the beginning" when introducing a correction or problem 
variation; it is only required that information which is actually chang- 
ing be respecified. Second, it should be noted that while data can be 
modified at any time, the effects of a modification on finite-difference 
equations and solution components can be realized only as a result of 
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explicit user requests. As indicated earlier, EPS does not actually 
assemble data to produce a solution until expressly directed to do so. 
This means that the user can make corrections as he sees fit during 
the course of a problem description without unnecessary loss of compu- 
ter time, but it also means that he must remember to be explicit after 
introducing a problem variation lest he be misled into thinking that 
he has a new solution when he has not. A misimpression of this sort 
is made possible by the fact that "current" solutions are deleted by 
EPS only when major problem revisions, namely, changes to the lattice 
structure in i„ "'-space or to the number of differential equations n, 
are requested. Thus, to consider a problem variation, introduce all 
changes of interest, then type form- and re 1 ax-commands , as described 
in Section 6, to obtain a new solution. If no major changes are intro- 
duced, the relaxation routine will use the previous solution as its 
starting point; it will also use previous values for parameters "omega' 
"delta" and "limit" unless new ones are specified. 

8.1. Modification of lattice structure in -£,j'-space. Data 
relating to lattice structure in i ,j -space can be modified in three 
ways : 

(1) New contours can be defined, e.g., to introduce holes in 

an existing domain, by typing additional append- and close- 
commands, as discussed in Section 3.4. 

(2) An existing contour curvename can be deleted in its entire- 
ty, e.g., to eliminate a hole or facilitate some other 
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gross lattice revision, by typing 

delete aurvename 
The computer's response to such a command has the form 

'CURVENAME' HAS BEEN DELETED. 

(3) An arc of contour aurvename, beginning at £,,,7',, passing 

through some intermediate point i ,j , and ending at t,,j,, 
can be deleted by typing 

delete i-,,0-,, £ 2 ,J 2 > ^ ,j' from aurvename 

Then a replacement arc can be specified and aurvename reclosed 
by typing additional append- and close-commands. For exam- 
ple, a new grid line, i=8 , can be introduced in a lattice 
bounded by a single contour "alpha" with vertices 0,0, 7,0, 
7,7 and 0,7 by means of the sequence 

delete 7,7, 7,2, 7,0 from alpha 
append 7,0, 8,0, 8,7 to alpha 
close alpha $ 

The computer's response in this case would be 

DELETION FROM 'ALPHA' HAS BEEN COMPLETED. 
ADDITION TO 'ALPHA' HAS BEEN COMPLETED. 
'ALPHA' HAS BEEN CLOSED. 

The user should note that 

fa) none of the i, , 1, of the delete-command need neces- 
sarily coincide with nodes which are contour vertices 
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in the -i, j -plane; 

(b) when the append -command is used to append a new arc 
onto an existing curve, it is necessary that the new 
arc and the existing curve have at least one end point 
in common; and 

(c) any boundary segment numbers associated with deleted 
arcs must be respecified after replacement arcs have 
been introduced. (Replacement arcs are automatically 
assigned a segment number of zero.) 

8.2. Parameter modification. The modification of quantities 
such as 

(1) lattice parameters x and y, 

(2) equation parameter n and the a(k,l) , b(fc,Z), .... g(k,l) 
and h (k) , 

(3) boundary condition parameters p(k,l,m), q(k,l,m) and r(k,m), 
and 

(4) relaxation parameters omega, delta and limit, 

as well as any auxiliary parameters which the user may have introduced, 
can be accomplished in a natural and obvious manner, namely, by reusing 
the set- and def ine-commands. For example, if it is decided that coef- 
ficient a (2,2) should assume the constant value 2.5, type 

set a(2,2)=2.5$ 
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regardless of how a(2,2) was specified in the past. The older defini- 
tion will be overwritten without any explicit computer response. Users 
with advanced applications or individuals wishing to use EPS as an on- 
line programming language should note the following points: 

(1) Information relating to parameter type , such as the number 
and range of subscripts (DIMENSION information, in the FOR- 
TRAN sense) or the distinction between constants and varia- 
bles, is updated automatically whenever a parameter is re- 
specified, and it is permissible not only to change a varia- 
ble into a constant or vice-versa, but also to change a sca- 
lar into an array or vice-versa. 

(2) Parameters cannot be deleted entirely in current versions 
of EPS, but significant amounts of storage can be recovered 
by using set-commands to change unneeded variable or sub- 
scripted parameters into scalar constants. Obviously, such 
steps are essential only in applications where storage limi- 
tations are a problem. 

8.3 . Modification of segment numbers. Like parameters, seg- 
ment numbers can be modified simply by reusing the command by which 
they are specified in the first place. Thus to assign a new segment 
number, type an appropriate impose-command, as described in Section 
5.2. The old segment number (or numbers) will be overwritten without 
explicit response. 
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9. Macros 

9.1. Introduction. In order to make possible the abbreviation 
of frequently used input sequences, EPS provides a macro facility simi- 
lar to those of many machine language assemblers and some high-level 
compilers. With EPS, a macro, or input substitution rule, can be speci- 
fied in terms of an identifier maoroname, a set of formal parameters 
{p^}, a set of separators {s^}, and a definition d(p-^;p2i - • • jPn^ as 
follows : 

let maoroname (p-^ s^ p 2 • •• s #_i P#) mean 
begin 

d(piiP2l- • -JPn^ 
end 

(The notation d(p-^jp 2 J ■ - ■ SPiy) is not meant to imply any particular input 
form, but merely represents an arbitrary input sequence in which the 
formal parameters (may) appear.) The computer responds to such a com- 
mand by printing 

INDICATED SUBSTITUTION FOR 'MAORONAME 1 WILL BE MADE. 

and, in processing subsequent input, replaces each call of the form 

maoroname (a-^ s^ a 2 ... s#_i <z#) 

where the au are actual parameters, or arguments, by an instance of the 
macro definition 



^This is the normal form of the let -command. Two minor varia- 
tions, discussed in sec. 9.4, are also permitted. 
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d(a 1 ;a 2 ; ■ ■ . ;a N ) 

in which each occurrence of each p^ is replaced by a literal copy of the 
corresponding a^. 

Before going into further details, we consider a specific example, 
namely, a macro to facilitate the specification of boundary values for 
problems with one dependent variable (see Section 5). Such a macro can 
be specified in terms of an auxiliary parameter "m" as follows: 

let prescr i be (va 1 ue along arc) mean 
beg i n 

set m=m+l , p(l,l,m)=0, q(l,l,m)=l 
define r(l,m)=value 
impose m along arc 
end 

Then, for example, typing a call such as 

prescribe(sqrt (fx) along 0,0, 2,0, 9,0) 

becomes completely equivalent to typing the command sequence 

set m=m+l, p(l,l,m)=0, q(l,l,m)=l 
define r (l ,m)=sqrt (fx) 
impose m along 0,0, 2,0, 9,0 

(It is necessary, of course, to assign "m" an initial value before 
making the first call.) Comparing this example with the general form 
of the preceeding paragraph, we see that the words "value" and "arc" 
are used as formal parameter names and the word "along" as a separa- 



-58- 



t 
tor. In general, formal parameters and separators may be arbitrary 

§ 
elements, other than parentheses and dollar signs. Thus, highly read- 
able constructions, such as the one given here, can be easily achieved. 

As the reader will recognize, macros are basically functions which 
operate in a somewhat broader domain than the functions of ordinary 
algebra. In fact, the only essential difference between macros, as 
implemented in EPS, and ordinary functions, stems from the fact that 
macro arguments are substituted literally, without binding or pre-evalu- 
ation. Thus, for example, the macro specified by 

let f(x,y) mean begin x/y end 

is not totally equivalent to the function given by 

f(x,y) = x/y, 

since a call for the macro made in terms of arguments "a+b" and "c+d" is 
interpreted according to 

f (a+b, c+d) •*■ a+b/c+d 

whereas the corresponding function value is 

f (a+b, c+d) = (a+b) I (c+d). 

Accordingly, certain precautions are necessary when dealing with macros 



The fact that "along" also appears in the macro definition has no 

significance. 
§ 
Words, numbers, and printed delimiters. See sec. 15.1. 

For other considerations, see sec. 9.5. 
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having algebraic significance. For example, to insure that macro "f" 
behaves like an ordinary function in all contexts, it is necessary to 
use the more explicit definition given by 

let f(x,y) mean begin ((x)/(y)) end 

In point of fact, it is generally recommended that macros not be used 
to represent single algebraic functions, but instead be used to repre- 
sent complete commands or command sequences, as in the macro "prescribe" 
shown earlier. In such cases, unnatural precautions are seldom neces- 
sary, for arguments tend to be complete expressions, lists of expres- 
sions, or whole commands, rather than single variables to be substi- 
tuted into larger expressions. 

9.2. When macros can be used. Like all EPS commands, the let- 
command can be employed whenever the user chooses, and thus macros can 
be specified as they are needed. Moreover, once specified, a macro can 
be called at any point, so long as the resulting substitution is appro- 
priate at that point. It is the user's responsibility to specify and 
call macros in such a way that well-formed and contextually valid sub- 
stitutions are produced. EPS, for its part, is blindly obedient to 
the user's instructions, and can easily be induced to make substitutions 
which are improper. In general, no real damage or computational loss 



Algebraic functions are best represented as auxiliary variables, 
specified through use of the def i ne-command. As currently implemented, 
this command does not admit to the explicit use of formal parameters, 
but, when used in conjunction with the set-command, does provide an effi- 
cient mechanism for evaluating formulae with varying parameters. 
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is incurred in such cases, but the user should understand that, since 
macro calls are literally replaced by their corresponding substitutions 
before command processing is initiated, error messages resulting from 
incorrect usage are normally expressed in terms of information which 
has been substituted for macro calls, rather than in terms of the calls 
themselves. Consequently, it is possible to receive error messages 
which, at the surface, appear totally unrelated to actual console input. 
A simple example is provided by the I/O sequence which follows: 

PROCEED: 

let pf(z) mean begin prnit z using bed -8f9.4*- end$ 

INDICATED SUBSTITUTION FOR 'PF' WILL BE MADE. 

PROCEED: 

pf( ux(2,5,7),uy(2,5,7))$ 

PRNIT IS NOT A COMMAND. TYPE 'LIST C0MMANDS$' FOR 

A REMINDER. 

Here the input line "pf( ... )$" has resulted in a somewhat oblique 
error message stemming from the earlier mistyping of the word "print". 
In this case the error message can be easily understood, but, as the 
adventurous or careless user will discover, more puzzling diagnostics 
are possible. 

9.3. Macro definitions. As the last example suggests, no formal 
checking of the definition portion of a let-command is performed until 
an actual call for the definition is introduced. Consequently, there 
are no explicit restrictions as to what a macro definition may contain, 
although it is implicit in the manner in which input is analyzed that 
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a definition may not contain a terminating "$", nor may it contain the 
sequence delimiters "begin" and "end" except in matched pairs. To illus- 
trate the varied forms which macro definitions may assume, it is useful 
to consider several specific possibilities. 

Definitions with missing formal 'parameters. It is permissible to 
introduce macro definitions in which some (or all] of the macro's formal 
parameters are unrepresented, as in 

let compute(x where y) mean begin print x end 

It is understood that in processing subsequent calls, actual arguments 
corresponding to the missing parameters are to be ignored. Thus, given 
the definition above, the call 

compute( m2+m3 where m2 and m3 are in ft. -lbs. ) 

is equivalent to 

print m2+m3 

It should be noted, in particular, that macro definitions may be total- 
ly degenerate. Thus the command 

let comment (x) mean begin end 

specifies a general mechanism for introducing explanatory remarks: 

comment(this entire call would be ignored) 
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Macros which call other macros. It is permissible to define a 

macro a which calls another macro 3 which calls another macro y, and 

+ 
so on, without restriction. For example, a user might type 

let pf(z) mean begin print z using bed -8f9.^ A " end 
let psc(t) mean begin pf (si n (t) ,cos (t) ) end $ 

Afterwards, he would find the call "psc(theta)" completely equivalent 
to the command 

print sin(theta) ,cos(theta) using bed -8f 9.4*- 

It should be noted that macro calls which appear within a macro defi- 
nition are not processed when the definition is initially specified. 
Consequently, the order of specification of a system of macros is im- 
material. Moreover, the effective meaning of a macro a which calls a 
macro B can be modified at any time by respecifying 3 (see Section 9.6). 

Macros with conditional phrases. It is possible to specify that 
an input phrase p be processed if and only if an expression e has a 
non-zero value by typing 

if e insert begin p end 

Such a sequence, called a conditional phrase or if -phrase, is permitted 



+ 
As a safeguard against infinite looping, EPS interrupts input 

processing whenever more than twenty macro calls are encountered and 

expanded without producing a single executable command. In such 

situations, the user is requested to supply explicit input indicating 

whether or not he wishes further substitutions to be performed. 



-63- 



anywhere, but can be used to particular advantage in a macro definition, 
as is illustrated by the following command: 

let zero(f ,df ,x,xo,tol) mean 
begin 
set x=xo 
print x,f 
if abs(f) grt tol insert 

begin 

zero(f,df,x,x-(f)/(df),tol) 

end 
end 

This command specifies a recursive macro which uses the tangent (or 

t 
Newton-Raphspn) method to find a iexo for an expression "f" with inde- 
pendent variable "x" and first derivative "df'J. Called with a start- 
ing value of "xo" for the independent variable, the macro produces a 
sequence of new values by calling itself repeatedly until the absolute 
value of "f" comes to be less than or equal to "tol". 

It should be noted that conditional phrases need not represent 
entire commands or command sequences, but may represent parts of com- 
mands as well. Consequently, the occurrence of an "if" in EPS input 
is not taken by the computer as a signal to terminate the command (if 
any) which preceeds it. Moreover, the computer must necessarily eval- 
uate the expression on which an if -phrase depends before it can initi- 
ate processing of the command which immediately preceeds that phrase. 
Therefore, in defining recursive macros it is necessary to avoid con- 



See, for example, Hildebrand [10] , sec. 10.8. 
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structions such as 

set k=k+2 

if k leq kmax insert begin phrase end 

where the command immediately preceeding the if-phrase affects the way 
the if-phrase is interpreted. This can always be done by inserting a 
stop-command immediately before the if-phrase, as in 

set k=k+2 

stop 

if k leq kmax insert begin phrase end 

The stop-command has no explicit computational purpose, but insures 
that the set-command is executed before the conditional phrase is en- 
countered. 

We conclude this section with a final example, namely, the command 
that was used to prespecify the EPS do-macro introduced in Section 7.1: 

let do(s for v=vl step dv until v2) mean 
beg i n 
stop 
if (dv)*((v2)-(vl)) geq insert 

begin 

set v=vl s 

do(s for v=v+dv step dv until v2) 

end 
end 

The reader should note again that this macro can be used to cause repe- 
tition of any command sequence "s" (including a sequence containing 
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other calls for the do-macro). He should also note that the parameters 
"vl", "dv", and "v2" represent arbitrary algebraic expressions and that 
the step size "dv" may be either positive or negative and either inte- 
gral or fractional. Less general versions of this macro, e.g., versions 
which assume a fixed step size or prohibit arbitrary expressions for 
"vl" and "v2", are obviously possible, and may be expected, in general, 
to function more efficiently. 

9.4. Special forms. The preceeding sections deal exclusively 
with macros specified by means of commands of the form 

let maoponame (p ^ s^ P2 ••• s ff-l Pff) mean 
begin 

end 

Actually, two variations of this form are permitted. One variation, 
used to specify parameter less macros , is given quite simply by 

let macroname mean begin d end 

This command causes the fixed definition d to be substituted for subse- 
quent occurrences of the identifier maovoname, and, while useful in many 
situations (see Section 9.7, Example 4), requires no detailed explana- 
tion. 

The second variation is used to specify self-repeating macros , and 
is given in terms of the constituents of a normal let-command, plus a 
"trailing separator" s^, as follows: 
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let maovoname (p± S\ P2 ••• S N-1 PN s N' mean 
begin 

end 

The presence o£ the trailing separator is taken to indicate that subse- 
quent calls may contain more than one set of arguments, each set being 
separated from the next by an Sjj, as in 

maovoname (a± s± a 2 — S N-1 a N S N 
b 1 s 1 b 2 ... s^ b N s N 



2 1 s l 2 2 " ' * S N-1 Z N^ 

The computer replaces each such call by a sequence of substitutions, one 
for each set of arguments, obtaining 

d(a-^;a 2 ; . . . ;a N ) d(b^;b 2 ; . . . ;b N ) ... d(z±;z 2 ; . . • ;zjf) 

A concrete example is provided by the command 

let pu(a,b,) mean 
begin 

set i=a, j=b print x, y, u(l,a,b) 
end 

This defines a macro which, for example, could be called with 

pu( 1,1, 2,J», 5,9) 

to cause execution of the command sequence 
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set 
set 
set 



=1, j=l print x, y , u(l,l,l) 
=2, j=4 print x, y, u(l,2,4) 
=5, j=9 print x, y, u(l,5,9) 



Note that no difficulty arises from the fact that the same element, 
a comma, is used for both intermediate and trailing separators. 

As the last example may suggest, self-repeating macros are extreme- 
ly useful in situations where repetition of a given command sequence is 
desired for sets of arguments which cannot be generated algorithmically . 
One such situation of general interest concerns the specification of 
numerical data in the form of an array. Suppose, for example, that 
the components of a vector "z" represent experimental values which must 
be specified manually and individually. Without macros such an opera- 
tion can be quite tedious, for it requires the explicit naming of each 
component, as in 

set z(l)=2.51, z(2)=3.7A, z(3)=4.97, z(4)=6.02, 
z(5)=7.76, z(6)=9.93, ... 

However, by typing 

let f i 1 1 z (x , ) mean begin set z(k)=x, k=k+l end 
one can eliminate the need to name each component, and write instead 

set k=i fillz(2.51, 3-74, 4.97, 6.02, 7-76, 9.93, ...) 

More generally, one can specify a macro for filling any vector "v" with 
a list of values "x", starting with component "v(ko)", by typing 
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let f i 1 1 vector (v from ko with x) mean 
begi n 

let f i 1 1 v (y , ) mean begin set v(k)=y, k=k+l end 
set k=ko f i 1 lv(x) 
end 

Similarly, one can specify a macro for filling an "n" column matrix 
"w" with values "x", assuming these are listed row-wise starting with 
component "w(l,l)", by typing 

let fi11matrix(w thru n with x) mean 
begin 

let fillm(y,) mean 
begin 
set if j grt n insert begin i=i+l, j=l, end 

w(i ,j)=y, j=j+l 
end 
set i=j=l fillm(x) 
end 

The reader should note that these last two macros have definitions 
which themselves specify macros. Nesting of macro definitions is, in 
fact, permitted to any depth. 

9.5. Separators, arguments, and imperfect calls. In several of 
the preceeding examples, argument separators have been chosen in such a 
way as to produce calls which read almost like English sentences. This 
is often a worthwhile practice, as it makes the order and significance 
of arguments easier to remember. However, it is important to recognize 
that the primary function of a separator is to provide the computer with 
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an unambiguous means of distinguishing one argument from the next. To 
do this, separators must be chosen with forethought, in accordance with 
the following restriction: a separator s, should not be any element 

which might occur in the preceeding argument a, , except possibly with- 

f 
in matching left and right parentheses. Thus, for example, if it is 

intended that a, be a list of algebraic expressions separated by commas, 
then s, should not itself be a comma. On the other hand, if a, is ex- 
pected to be a single algebraic expression, then s, could be a comma, 
since any commas which might occur in a, would necessarily be "protect- 
ed" by parentheses, as in "a(2,2)/fy". 

The protective nature of parentheses in macro arguments means 
that arguments themselves are subject to a restriction, namely, the 
requirement that they not contain unbalanced parentheses. In princi- 
ple, this is not a serious restriction, as it merely dictates that 
arguments be syntactically complete algebraic expressions, commands, 
macro calls, etc. In practice, however, difficulties can arise, for it 
is easy to introduce mismatched parentheses inadvertently when typing 
complicated expressions. The important point here is that mismatched 
parentheses cannot be detected as such by the EPS macro processor, for, 
at the level of the macro processor, calls with mismatched parentheses 
are indistinguishable from certain imperfect calls , which are admitted 
deliberately. As a result, mismatched parentheses simply lead to im- 
proper argument substitutions, which, in turn, generally yield ill-formed 



t 
Macro calls are analyzed by means of a left-to-right scan. There- 
fore, it is not strictly necessary, when selecting a given separator, to 
consider the possible forms of the argument which follows it. 
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or otherwise invalid commands. These are normally reported without 
serious consequences, although, as indicated in Section 9.2, they may 
be reported in terms which are somewhat obscure. 

Macro calls with argument lists which are incomplete by virtue 
of missing outer parentheses or missing arguments and separators are 
called imperfect calls, and are admitted by current versions of EPS 
in order to enable the user to take short cuts in certain restricted 
but common situations. Three cases that may be of general interest 
are discussed below. 

(1) It is permissible to omit the opening left parenthesis of 

a macro call's argument list provided (a) the first argument 
a± does not itself begin with a left parenthesis and (b) the 
call is not itself an argument to another call. Thus, for 
example, no difficulty is incurred by typing 

seq-^ do print z(k) for k=l step 1 until 8) seq2 

provided seq-, and seqj ^° not represent parts of another macro 
call in which the call "do print ..." is an argument. 

(2) It is permissible to omit the closing right parenthesis of 
an argument list provided this list is the very last thing 
in an input sequence and, therefore, is followed immediately 
by a "$". (In effect, the computer fills in as many right 
parentheses as are needed to balance left parentheses occur- 
ring earlier in the sequence.) This case and case (1) are 
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illustrated concurrently in the sequence which follows: 

set n=3, b=c=d=f=h=0 

do do(set a(k,l)=e(k,l)=g(k,l)=0 for 1=1 step 1 

until n) for k=l step 1 until n) 
do set a (k,k)=e(k,k)=g (k,k+l)=l for k=l step 1 

unt i 1 n $ 

Note, in particular, that it has been possible to omit both 
outer parentheses in the final call for macro "do". This 
sequence, incidentally, specifies governing equation parame- 
ters corresponding to the system { V u -. + u„ = 0, V u„ + 
u~ = 0, V 2 u^ = }, which is equivalent to the sixth-order 
equation V u, = (see Section 4) . 

(3) It is permissible, in general, to omit any or all of the argu- 
ments in a macro call and also to omit separators which would 
otherwise be left trailing at the right end of an argument 
list as a result of the omission of arguments. The computer 
uses formal parameter names in lieu of all missing arguments. 
Thus, by using parameter names which themselves have signifi- 
cance, one can specify macros with optional arguments. For 
example, a looping macro analogous to the FORTRAN DO-statement 
is specified by the command 

let dof(s for i=0,n,l) mean 
beg i n 

do(s for i=0 step 1 until n) 
end 
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Possible calls and corresponding substitutions for this macro 
are illustrated below: 




It must be emphasized that separators cannot be omitted when 
they occur at intermediate points in an argument list. For 
example, a call for macro "dof" which omits only the argument 
corresponding to parameter "n" must have the form 

dof (seq for u=u ,,5u) 

Similarly, a call which omits only the arguments correspond- 
ing to "i" and "n" must have the form 

dof (seq for=i 1 , ,&i) 

For maximum convenience, macros should be specified so that 
optional arguments occur in order of increasing likelihood 
of omission. 

The discussion of imperfect calls given here is not complete, but can 
be easily supplemented by direct experimentation. To facilitate exper- 
imentation, EPS provides a special output command, normally written in 
terms of a sequence p containing macro calls, as follows: 
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expand beg i n p end 

The computer responds to this command by printing its interpretation 
of the sequence p. For example, after the specification of macro "pre- 
scribe" shown in Section 9.1, the computer would respond to 

expand begin prescribe fx+fy along 5,5, 5,6) end $ 

by printing 

SET M=M+1, P(1,1,M)=0, d(l,l,M)=l DEFINE R(l,M)=FX 
+FY IMPOSE M ALONG 5,5, 5,6 

Note that it is not assumed that sequence p represents an executable 
command sequence, nor is any attempt made to execute p if it turns out 
to be an executable command sequence. Consequently, the expand-command 
cannot be employed to determine the significance of calls for recursive 
macros, such as the do-macro, whose proper operation hinges on actual 
command execution. Tt may, on the other hand, be employed to cause 
printing of arbitrary text strings and, therefore, is useful not only 
for testing (non-recursive) macros, but also for printing non-numeri- 
cal output, e.g., table headings and messages of confirmation (see Sec- 
tion 9.7, Example C) . 

9.6. Macro deletion and respecification. It is possible to 
delete one or more macros by typing a command of the form 



remove begin macroname , maoroname . 



end 
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For each identifier listed (and actually found to name a macro), the 
computer will respond with a message of the form 

MACRO 'MACRONAME' NOW INACTIVE. 

Subsequent uses of all identifiers will then be taken literally. Since 
the number of storage registers required to represent a macro is rough- 
ly twice the number of elements in the let-command used to specify the 
macro initially (the do-macro requires about 100 words of storage), it 
is worthwhile to delete macros, thereby recovering storage for other 
purposes, whenever they are no longer needed, 

A macro can be respecified at any time simply by introducing a 
new let-command. When this is done, the computer simulates an appro- 
priate remove- command before processing the new macro description. 
Therefore, it prints a "MACRO INACTIVE" message before printing the 
usual "INDICATED SUBSTITUTION ... WILL BE MADE." As may be evident, 
no more than one macro with a given identifier is permitted. Therefore, 
any let -command specifying a maoroname which is already in use neces- 
sarily forces the removal of the current description, even if the new 
description involves a different number of formal parameters or a dif- 
ferent set of separators. 

9.7. Additional examples. In this final section, three addi- 
tional examples, related specifically to the solution of boundary-value 
problems, are discussed. (See also Sections 10 and 11.) 

A. Specification of boundary conditions. The usefulness of macros 
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in specifying boundary conditions has already been indicated in terms 
of the macro "prescribe" shown in Section 9.1. This macro, however, 
is concerned with a very simple class of problems and boundary condi- 
tions. A more complicated situation and more imaginative use of macros 
is illustrated by the command sequence below, which specifies a system 
of macros for use in the solution of elasticity problems (see Section 
4.3, Example D, and Section 5.5). 

let make (a=b ,c=d along e) mean 

beg i n 

set m=m+l, p (l , 2 ,m)=p (2 ,1 ,m)=q (l , 2 ,m)=q (2 ,1 ,m)=0 , 
a=b, c=d 

impose m along e 

end 
let dx mean begin p(l,l,m)=0, q(l,l,m)=l, r(l,m) end 
let dy mean begin p(2,2,m)=0, q(2,2,m)=l, r(2,m) end 
let tx mean begin p(l,l,m)=l, q(l,l,m)=0, r(l,m) end 
let ty mean begin p(2,2,m)=l, q(2,2,m)=0, r(2,m) end 

This sequence enables the user to specify boundary conditions in terms 
of the parameterless macros "dx", "dy", "tx", and "ty", which corres- 
pond, respectively, to x-displacement , ^-displacement, x-traction, and 
y-traction conditions. For example, to indicate that the boundary arc 
5,-2, 6,-2, 23,-2 is to have no x-displacement and is to be subjected 
to a y-traction of one unit, the user could type 

make(dx=0, ty=l along 5,-2, 6,-2, 23,-2) 

Obviously it is possible to formulate similar systems of macros to facil- 
itate the specification of boundary conditions in any application area. 
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B. Numerical differentiation. As indicated in Section 7, EPS 
provides built-in functions "ux", "uyy", etc. for calculating approx- 
imate values for first and second derivatives du,/dx, 9 u^/dy , etc. 
based on current states of solution matrix "u" and lattice parameters 
"x" and "y". At times, however, the user may find it necessary to 
supplement these built-in functions with procedures of his own which, 
for example, calculate derivatives of the u^ with respect to variables 
other than x and y, or which calculate derivatives of quantities other 
than the u, . Here a procedure is considered which calculates approxi- 
mate slopes for an arbitrary function / in the directions of the grid 
lines of a lattice. It is assumed that values of f for adjacent nodes 
of the lattice are well-behaved with respect to changes in arc length 
s measured along the grid lines, justifying the representation of f(s) 

in the neighborhood of any three consecutive nodes with a second-order 

+ 
polynomial. By choosing coefficients so that the polynomial matches 

given values /. [k = 1,2,3] for arguments s< corresponding to the three 

nodes and then differentiating, one obtains the following approximation 

for the slope of / at intermediate node s~' 

\df\ (s 2 - s 1 ) 2 (f 3 - f 2 ) + (s 3 - s 2 )V 2 - f ± ) 

3sj _ (s 2 - s 1 )(s 3 - s 2 )(s 3 - s 1 ) 

This approximation has been incorporated in the following command se- 
quence, which specifies a macro that causes printing of the approxi- 



In general this assumption demands, among other things, that 
the straight line segments which constitute the grid lines change di- 
rection slowly, i.e., the grid lines must be reasonable approximations 
to smooth curves . 
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mate slope in the '^"-direction, where "i" or "j" must be substituted 
for "k", of any quantity "f" at any intermediate node "i2,j2": 

let slope(f ,k, i2, j2) mean 
beg i n 
set i=i2, j=j2, f2=f, x2=x, y2=y, 

k-k-1, fl=f, ds21=dsfcn, 

k=k+2, f3=f, ds32-dsfcn 
print dfds 
end 
define dsfcn=sqrt ( (x-x2)power 2+(y-y2)power 2), 
dfds=(ds21*ds21*(f3-f2)+ds32*ds32*(f2-fl))/ 

(ds21*ds32*(ds21+ds32)) 

As a result of this sequence, the approximate slope of solution compo- 
nent u 2 in the ^-direction at node (5,6) could be obtained by typing 

s1ope(u(2 f i,6),i,5,6) $ 

Note that the first argument in this call could also be written 
"u(2,i,j)'\ but not "u (2,5,6)", the latter being constant for all "i". 
To avoid confusion and minimize typing when several calls are antici- 
pated, it is advisable to define as auxiliary parameters functions to 
be used as first arguments. For example, type 

def ine psi=u(2,i ,j) slope(psi , i ,5,6) $ 

All subsequent calls for derivatives of w- could then be made in terms 
of "psi". 

C. Numerical integration. In many problems, integral properties 
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of a solution or its derivatives are of interest. EPS provides no built- 
in mechanism for computing integrals, so the user must supply his own. 
For most purposes, integrals can be adequately approximated by use of 
the trapezoidal rule, which is formulated in terms of an integrand f(z) 
and a set of arbitrarily spaced sample points z^, s m+1 » •••» z n as fol ~ 
lows : 



z n n 



$ f{z)dz * I \ [f(z t ) + f (3 £ _ x ) ] (a £ - »£_]_) 



m 



Shown below is a command sequence implementing this formulation for 
integrations along grid lines of constant "k" between nodes "k=ko, 
l=m" and "k=ko,l=n". 

def i ne news i g=s i g+ . 5* (f 1+f 0) * (zl-zO) 
let sum(f,z for k=ko, l=m,n) mean 

begi n 

set k=ko, l=m, f0=f, z0=z, sig=0, dl=l sign(n-m) 

do( set fl=f, zl=z, sig=newsig, fO-fl, z0=zl 
for l=m+l step dl until n) 

expand begin integral of f wrt z along k=ko, 
l=m,n is end 

print sig 

end 

With macro "sum" one could, for example, obtain the approximate value 
of / (dujdx)dy for a path C corresponding to the grid line i=5 , for 
j=0 through 20, by typing 

sum(ux(l,5,j),y for i=5, j=0,20) $ 
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This call would result in a computer response of the form 

INTEGRAL OF UX(l,5,J) WRT Y ALONG 1=5, J=0,20 IS 

±X.XXXXXXXE±XX 

As in Example B, it is advisable to use auxiliary parameters for argu- 
ment "f" if multiple calls are anticipated. Various generalizations 
(e.g., macros to perform area integrations) and refinements (more accu- 
rate integration procedures) are obviously possible. 
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10. Iterative problem solving 



10.1. Introduction. This section discusses iterative techniques 
for solving boundary-value problems for nonlinear partial differential 
equations, for treating free-surface problems, and for constructing po- 
tential flow nets. In all cases 



the basic strategy is as illus- 
trated in Figure 9. First, a 
finite-difference lattice and 
differential system approximat- 
ing the situation of interest 
are specified and an initial so- 
lution is obtained by means of 
the procedures described in Sec- 
tions 3 through 6. Then equa- 
tion or lattice parameters are 
revised on the basis of the ini- 
tial solution to establish a 



(begin^) 



SPECIFY INITIAL 
PROBLEM 



REVISE 
DATA 



OBTAIN NUMER- 
ICAL SOLUTION 



fail 




Fig. 9. Iterative strategy for 

nonlinear equations, etc. 



more accurate representation of the true situation, and a second solu- 
tion is obtained. This process is repeated until convergence as defined 
by some suitable criterion is achieved. 

Detailed procedures for several specific cases are given in Sec- 
tions 10.2 and 10.3. However, it is emphasized that the basic strategy 
of Figure 9 is applicable to many situations not even implicitly consi- 
dered. For example, engineering design problems requiring the optimi- 
zation of part shape or other physical properties are obviously amena- 



ble to the same approach. At the same time, however, the procedures 
actually considered may not be suitable for every possible application. 
Convergence may be prohibitively slow in some cases, while outright 
divergence may be unavoidable in others. (Techniques which can some- 
times be used to accelerate convergence or increase numerical stabil- 
ity are discussed in Section 10.4.) In general, nevertheless, these 
procedures are remarkably dependable, and provide what is in many in- 
stances the only available means for realizing a solution. 

10.2. Nonlinear equations. The form of governing system (1), 
Section 4.1, coincides with that of various nonlinear systems when 
coefficients a, „, b,„, ..., g,„ are considered to be functions not 
only of the independent variables x and y , but also of the dependent 

variables u n , u», ..., u and their derivatives. The EPS difference 

1 2 n 

equation generator does not explicitly admit such coefficients, but 
it has been found that solutions can frequently be obtained by itera- 
tion in cases where no derivatives of order greater than one are in- 
volved, i.e., where 



. 8u- du 3u_, du , 

a ki = v-ki^^y > u i> • • • > u n > !.••■' n> !»•••> yj > 

dx dx dy dy 

, , , du-, du du-. du . 

kl = b kl^ X ' y,U l' ' ' ' ' U n' -'•'•' -' !'•••» yJ > 

dx dx dy dy 

etc. 



What this entails, basically, is the solution of a sequence of linear 
problems, each of which has spatially varying coefficients evaluated 



by using the previous solution of the sequence. To be specific, esti- 

o 
k 



mates u° = uUx,y) [k = l,2,...,n] for the behavior of the dependent 



variables are used to formulate an initial problem with coefficients 

-1 ^0^0 ~v, 

a kl = a kl (x ^' U l'--^ U n^'---'^' T J -'---' V Jl) ~ a kl ix > y) ' 

dx dx dy dy 

etc. 



The solution components u) = ul(x,y) obtained by solving this initial 
problem are then used to formulate a second problem with coefficients 

. 1 i du., 9w du n du . i , -. 

a kl = «fe£^.J/."i.---." n ._i.-".-_«._i.---.— S> = ^x t y), 

dx dx dy dy 

etc. , 



and so on, until finally the difference between one solution and the 
next falls below some acceptable tolerance. Convergence is not always 
guaranteed, of course, but in many practical applications it in fact 
proceeds at such a rate that only a few iterations are necessary. 

In terms of computer input, it is generally convenient to use 
the following procedure: 

(1) Specify coefficients a,„, b^n, etc. in terms of auxiliary 
parameters representing their various arguments, e.g., if 
a = a 2 Ax,u du /dy) = 8x + u^ - \du^/dy\, type 

define a (2 ,l)=8*f x+f ul-abs (f uyl) 

(2) Specify auxiliary parameters introduced in step (1) using 
the fit-function discussed in Section 4.2, making certain 
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that reasonable initial estimates for all pertinent depen- 
dent variables and derivatives are included. For example, 
if it is decided that w° = xy , hence du°/dy = x, type t 

define fx=fit(x), fy=fit(y). 
ful=fx»fy, fuyl=fx 

(3) After obtaining a first solution, redefine the auxiliary 
parameters that represent dependent variables in terms of 
appropriate components of the solution matrix, e.g., 

define ful=f it (u(l, i ,j) ) 

Similarly, redefine parameters representing derivatives of 
the dependent variables, in this case introducing auxiliary 
matrices into which values for the derivatives produced by 
means of the built-in functions "ux" and "uy" can be stored 
prior to each succeeding iteration, e.g., 

define fuyl=f i t (dyl ( i ,j)) 

[Coefficient parameters cannot be defined directly in terms 
of functions "ux" and "uy" because these functions cannot 
be referenced during the actual generation of finite-differ- 
ence equations.) 

(4) Specify pertinent components of auxiliary matrices intro- 



TThis example demonstrates the general case, but in practice it 
is normally convenient and numerically just as satisfactory to use con 
stant values for the u^ in formulating an initial problem. Generally 
it is appropriate to begin with some limiting case, e.g., ui = 0, cor- 
responding to a standard linearization for the problem of interest. 
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duced in step (3), if any, e.g., by typing input sequences 
such as 

do( do( set dyl(i,j)=uy(l,i,j) 

for i=0 step 1 until imax) 
for j=0 step 1 until jmax) 

(5) Type form- and re lax-commands to obtain a new solution. 

(6) Repeat steps (4) and (5) as needed until satisfactory con- 
vergence is achieved. 

Note that the typing involved in repeating steps (4) and C5), as well 
as in carrying out any computations needed to test for convergence, can 
be reduced considerably by means of appropriate macros. 

To illustrate the procedure outlined above, three specific exam- 
ples are considered. 

A. Heat conduction in a material with temperature-dependent 
conductivity. The steady-state temperature 9 in a material with con- 
ductivity fe(9) = (9/6 ) S feo, where 9 , s, and k are constants, satis- 
fies the following nonlinear equation: 



t 
It may also be desirable at this point to revise boundary con- 
dition or other parameters. Normally the effects of nonlinearities 
show up gradually as the inhomogeneities or forcing functions in a 
differential system are increased from zero. Thus numerical insta- 
bilities can sometimes be avoided if the quantities hfa in governing 
system (1) and the quantities rf^ in boundary equation (9) are ini- 
tially assigned relatively small values, then gradually increased as 
the iteration progresses. Similarly, multipliers of nonlinear terms 
in governing equations can sometimes be varied in order to control 
instability. See Example C. 
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div [(6/8o) S /c grad 0] - 0. + 



When expressed in two-dimensional form and in terms of Cartesian coor- 
dinates, this equation corresponds to governing system (1) with 

"1 -* e, 

n - 1, 
a ll = e ll = ("i/9o) S feo» 
fc ll " G ll " d\l - fn - 0ii * *! - 0. 

Accordingly, governing equation parameters are specified by means of 
a command sequence such as the following: 

set n=l, b=c=d=f=g=h=0 
define a=e*fe *(t/6 ) power e 

Auxiliary parameter "t" is then given some appropriate initial defi- 
nition and a first solution for the desired lattice and boundary con- 
ditions is obtained. After this, "t" is redefined in terms of solu- 
tion matrix "u", i.e., the command 

define t-f !t(u(l,i J)) 

is given, and form- and re lax-commands are issued repeatedly until 
suitable convergence is achieved. Convergence rates of about one 



This equation can be linearized (converted to Laplace's equa- 
tion) by means of a variable transformation. However, since linear 
boundary conditions relating heat flux and temperature become nonlin- 
ear when expressed in terms of the transformed variable (see Ames [11], 
pp. 21-23), it is often just as convenient to leave the equation in its 
original form. In any event, the procedure considered here is obvious- 
ly applicable in other situations where linearization by variable trans- 
formation is not possible. 
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decimal digit per iteration have been observed in several test cases 
for various values for exponent s between -1 and 3. 

B. Isentvopic compressible flow. The continuity equation for 
the steady, irrotational , isentropic flow of a perfect gas can be 

written in terms of a velocity potential $ and specific heat ratio 

t 
k as follows : 



div 



(l-|grad 0| 2 ) 1/(fe 1) grad $ 



0. 



When Cartesian coordinates are employed, the two-dimensional form of 
this equation corresponds to governing system (1) with 

"1 + $> 
n = 1, 

ft i a i2a/(fe-i) 

a ll = e ll = ( 1_ |g rad u ±\ ■> > 

h \\ = c ll = d \\ = f\\ = ^11 = h l = °- 

Thus equation parameters are specified by means of a sequence such as 

set n=l, b=c=d=f=g=h=0, kl=l/ (k-l) 
define a=e= (l-vv) power kl 

where auxiliary parameter "vv", corresponding to the quantity | grad $ j . 
is normally set to zero initially in order to obtain a first solution 
representing ordinary potential flow. Afterwards it is convenient to 



t 
This equation is obtained by noting that density p and veloc- 
ity v for isentropic flow are related by p/p = (l _ |v/v m axl ) ^ -1 ) , 
where p Q is the stagnation density and v ma x the velocity magnitude at 
zero enthalpy, and by defining potential $ so that grad $ = v/v m ax- 
See Shapiro [12], Arts. 4.3 and 4.4. 
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redefine "vv" in terms of a single matrix, as in 

define vv=f i t (vvm(i ,j)) 
where the components of "vvm" are to be computed in accordance with 

vvm(i.j) <- f|grad $| 2 = (3$/3a:) 2 +0*/3z/) 2 

[ Ut (f,j) 

(Note that this is more efficient than having two matrices, one for 
3$/3ar and one for 3$/3i/, though it represents a slight departure from 
the procedure defined earlier.) Thus, for example, with a lattice 
which is rectangular in ijj-space, say with vertices (0,0), (20,0), 
(20,15), and (0,15), a second solution is obtained by typing 

let revise mean 
begin 

do( do( set vvm(i , j )=»vvfcn 
for i=0 step 1 until 20) 
for j-0 step 1 until 15) 
end 
define vvfcn=ux(l, i , j)power 2+uy (l, i , j)power 2 
revise form relax$ 

After this, succeeding solutions are obtained simply by typing "revise 
form relax$" repeatedly, or, better yet, by introducing and calling 
a recursive macro which automatically causes this sequence to be re- 
processed until some appropriate convergence condition is satisfied. 
For example, the macro defined below causes iterations to be repeated 
until fractional changes in "vvm(icri t, jcri t)" fall below "tol", where 
presumably "(icri t , jcri t) M is a node of particular concern. 
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let repeat (tol) mean 
begi n 

set vvc=vvm(icri t,jcri t) 
revise stop 

if a bs ( ( vvm (icrit,jcrit)-vvc) /vvm (icrit,jcrit)) 
geq tol insert 

begin 

form relax repeat (tol) 

end 
end 

Convergence rates of one decimal digit every two or three iterations 
have been observed in several test cases where boundary conditions were 
chosen so as to produce Mach numbers of unity at one or two points in 
the flow. Interestingly, there appear to be no problems with numerical 
instability even when Mach numbers exceed unity by small amounts, though 
solutions may not be physically tenable in such cases. 

C. Visaous incompressible flow. The dynamical equations for 
steady, two-dimensional, viscous, incompressible flow can be written 
in terms of dimensionless Cartesian coordinates x and y, a dimension- 
less scalar vorticity function cj, dimensionless stream function ¥, and 
Reynolds number R as follows: 



9x 2 + 9z/ 2 *"l 2y dx dx dy 



2u+ + = ° 



t 



See, for example, Schlichting [13], pp. 58f. 
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It has been found that solutions to these equations for flows around 
obstacles can be obtained for modest Reynolds numbers, say, values 
less than 500, based on the characteristic dimension of an obstacle. 
However, to avoid numerical instabilities, it is generally necessary 
to proceed in stages, starting with the limiting case R - and grad- 
ually introducing larger values, e.g., R = 50, 100, 150, ..., until 
a solution for the value of interest is obtained. Note, first, that 
the nonlinear terms of the first equation can be rewritten in the form 

Thus a correspondence with governing system (1) can- be established in 

which 

u± -*■ w, u 2 -*■ V, 
n - 2, 

a ll = a 22 - e ll = e 22 " 1 » 021 " 2 ' 
b 12 = -i? «]_, i]_2 * R "1* 

and in which all other coefficients are zero. Equation parameters 
corresponding to this formulation are specified by typing 

set n-2, a(l,2)-a(2,l)-b(l,l)-b(2,l)-b(2,2)-c- 
d(l,l)-d(2,l)=d(2,2)-e(l,2)-e(2,l)-f-g(l,l)- 
g(l,2)-g(2,2)-h-0, a(l,l)-a(2,2)«e(l,l)« 
e(2,2)-l, g(2,l)-2 

define b(l,2) — rn*w, d(l,2)»rn*w 

where auxiliary parameters "rn" and "w", corresponding to R and u, 
respectively, are ordinarily set to zero before difference equations 
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are requested for the first time. After an initial solution is ob- 
tained, "w" is redefined by typing 

define w=f i t (u (l , i , j ) ) 

and succeeding iterations are performed by assigning "rn" a small non- 
zero value, typing form- and rel ax-commands until partial convergence 
is achieved, reassigning "rn" a slightly larger value, typing more 
form- and rel ax-commands , and so on. It has been observed that as 
larger values for "rn" are introduced, iterations display a growing 
tendency to wander and diverge. To counteract this tendency, smaller 
values for relaxation factor "omega" and techniques described in Sec- 
tion 10.4 can be introduced. Eventually, however, a Reynolds number 
is reached for which it is impossible to maintain numerical stability 
and an acceptable convergence rate simultaneously. 

10.3. Iterative lattice revision. In the analysis of seepage 
through porous media, elasto-plastic material behavior, fluid jets, 
and other phenomena, it is common to encounter boundary-value problems 
with free surfaces, i.e., with boundary segments which are initially 
undefined. Generally it is possible to associate with a free surface 
more than one set of boundary conditions. Hence the shape of the sur- 
face can usually be determined by first assuming some initial configu- 
ration and obtaining a solution satisfying one set of boundary condi- 
tions, then using the discrepancy between that solution and the remain- 
ing set of boundary conditions to revise the configuration, and so on. 

With EPS it is easy to carry out such a process if changes in 
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boundary geometry from one iteration to the next can be made without 
modifying lattice structure in i, ,j -space. Depending on the circum- 
stances, this can usually be accomplished by means of one or both of 
the following techniques: 

(1) The mapping between x,y- and i,j -coordinates can be expressed 
in terms of interpolation formulae giving the locations of 
interior nodes as functions of the locations of boundary nodes 
(see Section 3.5). Then iterations can be performed simply by 
revising auxiliary parameters which specify the locations of 
nodes lying on the free surface. Interpolation formulae obvi- 
ously must be devised so that grid lines cannot assume illegal 
configurations as the free-surface nodes are repositioned. 

(2) The mapping between x,y- and i ,j'-coordinates can be given in 
matrix form and each node adjusted independently to establish 
a lattice whose nodes lie along characteristic function con- 
tours. For example, in potential flow problems nodes can be 
made to lie along lines of constant potential and constant 
stream function, so that the lattice becomes a flow net. 

The second technique is sometimes useful even when no free surface is 
present, since the lattice it produces may have intrinsic utility and 
may also accord more accurate solutions. (Experiments indicate that 
this is generally true in the case of potential flow problems.) 

Because the strategy for revising lattice parameters varies from 
application to application, general input procedures corresponding to 
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the above techniques cannot be presented. Instead, procedures appro- 
priate to a specific but representative problem are described. Consi- 
der, in particular, steady seepage through a porous dam with trapezoi- 
dal cross section and horizontal, impervious base (see Figure 10). In 
this situation, the free surface is the interface between the region 
of seepage and the dry portion of the dam (curve AB in the figure) . A 



Free surface 
(approximate) 




Impervious base 



Fig. 10. Schematic of trapezoidal dam. 

boundary-value problem for the region of seepage can be formulated in 
terms of the piezometric head h, which under suitable assumptions plays 
the same role in seepage problems as the velocity potential in potential 
flow problems. In particular, if it is assumed that Darcy's law is sat- 
isfied and, further, that the permeability of the dam is the same in all 
directions and at all points, then one arrives at the following formula- 
tion (coordinates x and y , local normal ft, and dimensions h and h ? are 
as shown in Figure 10) : 



t 



See, for example, Polubarinova-Kochina [14], pp. 31-35. 
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Governing equation: —^ + — 2 = 0. (13) 

Boundary condition along 

dh 
free surface AB: t— = , h = v ; 



(14) 



dn 
exposed discharge face BC : h = y ; 

submerged discharge face CD: h = hj ; 

dh 
impervious base OD: -^— = 0: 

dn 

entrance face OA: h = h-i . 



For each of the above expressions, an equivalent expression can be given 
in terms of a stream function ¥, related to h by the Cauchy-Riemann 
equations 

ii = _ Ml and $1 = $h 

dx dy dy dx' 

Thus it is possible to arrive at the following alternate formulation: 



3 2 y 3 2 ¥ 

Governing equation: -r— 2 + -r— 2 = 0. (15) 

Boundary condition along 

free surface AB: y = v» = constant, ^- = &; 

A dn 34 



11 _ 

dn 


34' 


3¥ 
dn 


= 0; 


¥ = 


0; 


_3¥ 
3n 


= 0. 



(16) 



exposed discharge face BC: 
submerged discharge face CD: 
impervious base OD: 
entrance face OA : 

Here i> represents arc length along the boundary in the clockwise direc- 
tion. Note that along the free surface, which is fixed only at entrance 
point A, two boundary conditions are stipulated in each formulation. 
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Input procedures corresponding to the two iterative techniques 
outlined earlier and appropriate to the seepage problem just described 
are now considered. The first procedure, based on the use of interpo- 
lation formulae and aimed solely at the determination of the free sur- 
face and the head distribution h = h(x,y), is best suited to dams with 
steep discharge faces. The second procedure, where all nodes are re- 
vised independently to produce a flow net (hence solutions for both h 
and f) , is applicable to a wider range of configurations. Both proce- 
dures are outgrowths of an investigation carried out in collaboration 
with Lo [15], whose thesis may be consulted for further details and 
specific examples. 



A. Solution based on use of interpolation formulae. For dams 
with steep discharge faces it has been found to be convenient and appro- 
priate to use lattices of the type shown in Figure 11. Here each grid 




x ,1 



Lattice for seepage region of dam with 
interpolated nodal distribution. 
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line of constant i coincides with a straight line connecting a point on 
the base of the dam with the intersection point F formed by extending 
the entrance and discharge faces. For each such line, nodes are dis- 
tributed between a node (i,0) on the base and a node (^»J max ) on the 
free surface in a preestablished manner, so that as the free-surface 
node is moved up and down, the other nodes move by proportional amounts. 

iterations are initiated by estimating the elevation y~ i of 

J b u >^max 

each free-surface node and obtaining a numerical solution for h = h(x,y) , 
based on Eqs . (13) and (14), which satisfies the boundary condition 
dh/dn = along the free surface. Differences between this numerical 
solution and the estimated free-surface elevations--i . e. , discrepancies 
between the solution and the boundary condition h = y--are then used to 
compute new elevations, etc. It has been found that for a free-surface 
node where the discrepancy is 



hh . . = h. . - y • ■ , 
l »^max 1- '^max l '^max 



a new elevation y '. . can be computed using 

^ "?max 

y'. . = y . . + Ah. . sin 8^ , (17) 

l "?max ' 2, "'max l ' l 'max 

where Q^ is the angle which grid line i makes with the base of the dam 
(see Figure 11). When this expression was employed in various tests, 
the yj ■; converged at rates of one decimal digit every four or five 

iterations . 

In order to carry out the above procedure on the computer, the 
user can proceed as follows: 
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(1) Specify lattice structure in i,j -space in the usual manner 
(see Section 3.4], and define a mapping algorithm in terms of 
coordinates "xf" and "yf" of intersection point F, length "xd" 
of the base of the dam, vector components "yfree(i)" corre- 
sponding to the yi j , and interpolation functions "fi" and 
"fj" as follows: 

define x=xbase+(xf-xbase) -y/yf , xbase=f i -xd , 
y=f j -yf ree(i ) 

Note that it is permissible to express "x" in terms of "y" or 
vice versa, but not to use implicit formulations in which "x" 
and "y" both reference one another. 

(2) Specify parameters "xf", "yf", and "xd" in accordance with the 
known (fixed) geometrical features of the dam, and define "fi" 
and "fj" so as to produce the desired nodal distribution. As 
results are generally more accurate when nodes are crowded 
near exit point B of the free surface, it may be appropriate 
to type 

define f i= ( i /i max ) power pi, f j = (j/j' max ) power pj 

where pi and pj are positive numbers less than unity. (In 
Figure 11, pi - pj - 0.65.) 

(3) Specify initial estimates for the elevations "yfree(i)" of the 
free-surface nodes. (In many cases it is satisfactory to use 
values corresponding to a straight line connecting the known 
entrance point A to some assumed exit point B .) 
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(4) Specify the governing equation and boundary conditions shown 
in Eqs . (13) and (14), using, in particular, the condition 
dh/dn = along the free surface, by typing 

set a=e=n=l, b=c=d=f=g=h=0 
impose 1 along 0,0, 0,1, 0,j' max 
impose 2 along i max ,0 , i max ,l, W'^max 
impose at i max J max 
set p(l,l,0)=q(l,l,l)=q(l,l,2)=l, 
q(l,l,0)=p(l,l,l)=p(l,l,2)=0 
define r(l,2)=rr(fy grt hy) , rr (l)=fy=f i t (y) 
set rr(0)=h 2 , r(l,0)=0, r(l,l)=h 1 

Note that no impose -commands are given for the boundary arcs 
corresponding to base OD and free surface AB , hence the seg- 
ment number is assigned automatically by the computer. How- 
ever, it is necessary to assign segment number explicitly at 
node (^ maX jJ max ) , i.e., exit point B, to insure that the con- 
dition dh/dn = 0, not the discharge face condition h = y, is 
used. The computer elects "boundary value conditions" over 
"boundary slope conditions" when there is a choice, unless ex- 
plicitly directed to do otherwise. 

(5) Obtain an initial solution by using the form- and relax-com- 
mands as indicated in Section 6. 

(6) Specify a function for calculating new free-surface elevations 
in terms of parameters "s(i)" representing the sin 0^ of Eq. 
(17) by typing 

define newy=yf ree(i ) + (u(l, i ,j' max ) -yf ree ( i ))»s (i ) 
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Then specify the "s(i)" by typing 

do(set s(i )=yf/sqrt (yf power 2+ (xf -xbase) power 2) 
for i=l step 1 until i max ) 

(7) Revise lattice geometry by typing 

do(set yf ree ( i )=newy for i=l step 1 until i max ) 
and obtain a new solution by typing form- and re 1 ax-commands . 

(8) Repeat step (7) until adequate convergence is obtained. 

Note that the macro specified below could be called in lieu of steps 
(7) and (8) in order to cause iterations to be repeated until elevation 
changes at exit point B fall below "dy": 

let adjust (dy) mean 
beg i n 

set i=t max stop 
if abs (newy-yf ree ( i ) ) geq dy insert 

beg i n 

do(set yf ree(i)=newy for i =1 step 1 until i max ) 

form relax adjust(dy) 

end 
end 

Generally convergence is slowest at exit point B, and, moreover, eleva- 
tion errors are usually greatest in the vicinity of point B once conver- 
gence has been attained. However, surprisingly accurate results, e.g., 
maximum elevation errors of less than three percent of upstream water 
level h lt can be obtained with fewer than 100 nodes if the nodes are 
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appropriately distributed and if the discharge face is relatively steep, 

t 



say with 90° < 7 - < 135°. 
max 



B. Solution based on the generation of a flow net. In the proce- 
dure described above, the only objectives were the determination of the 
free surface and the head distribution h = h(x,y). Considered here is 
a procedure which generates an entire flow net, i.e., it produces a lat- 
tice whose grid lines approximate contours of constant h and constant 
stream function 'f. Note that the free surface itself is a contour of 
constant ^ (see Eq . (16)). An interesting feature of the present proce- 
dure is that the free surface is determined automatically, that is, the 
free-surface nodes are adjusted in exactly the same manner as nodes ly- 
ing in the interior of the lattice. Further, since the true contours of 
constant h and ^ are orthogonal, quadrilateral lattice elements are pro- 
duced which are virtually rectangular, a fact which is partially respon- 
sible for this procedure's increased accuracy. 

In brief, the procedure is as follows. An initial lattice having 
the topology of the desired flow net is specified and numerical solu- 
tions for h and ¥ are obtained which satisfy the boundary conditions h = 
y and 3¥/3n = dy/di, along the estimated free surface. (For convenience, 
h and V are treated as distinct dependent variables and solutions for 
both quantities, based on Eqs . (13)- (16), are obtained simultaneously. 
Obviously, either quantity could be determined from the other by inte- 
grating the Cauchy-Riemann equations; however, with EPS this requires 



"''See Lo [15] , p. 4-33. 
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more input and does not necessarily save computer time.) Since the ini- 
tial lattice will presumably not have the proper geometry, neither the 
additional free-surface conditions dh/dn = and Y = constant nor the 
requirement that (other) grid lines be lines of constant h and X V will be 
satisfied. Instead, there will be differences kh and A^ at most nodes 
between the desired and realized values of h and ¥. These differences 
and numerical values for the derivatives dh/dx , dV/dy, etc. are used to 
revise lattice geometry, then new solutions for h and V are obtained, 
and so on, until all requirements are satisfied within tolerable limits. 

A workable technique for revising lattice geometry involves the 
calculation of pairs of full corrections Ax and Ay for all nodes, then 
the replacement of the nodes' current coordinates (x,y) by new coordi- 
nates (x',y') given by 

x' = x + yAx, y' = y + dAy. (18) 

Here u is a relaxation factor which for best results is typically given 
a positive value less than unity at the outset, then gradually increased 
to unity as iterations progress. Full corrections Ax and Ay are deter- 
mined in the following manner: 

(1) For interior and free-surface nodes, where x and y can be ad- 
justed independently, the truncated Taylor expansions 

., dh . , dh . 
M = to ** + 3y" ^ 

A. = f Ax + f Ay =-fLAr + !£ Ay 

dx 6y y dy dx 

are inverted to give 
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Ax = — 



Ay 



1 

J 



dh .-, dh .... 

-r— Ah - t— A¥ 

dx dy 

dh . , dh .„, 

r— Ah + — AH 7 

dy dx 



(19) 



where 



9x 


2 

+ 


'8/zl 

, 3 .vJ 



(The Cauchy-Riemann conditions are invoked in the expansion 
for A^ to minimize the number of derivative calculations.) 

(2) For boundary nodes constrained to remain on prescribed arcs, 
say where x and y are given parametrically in terms of some 
variable t, i.e., where 

x = x(t) , y = y (t) , 

either Ah or A^ is implicitly taken to be zero and full cor- 
rections are calculated using 

Ax = x(t+At) - x(t), Ay = y(t+At) - y (t) , 



where 



At 



At 



dh dx dh_ dy 

dx dt dy dt 

dx dt dy dt 



Ah if Af = 0, 



-1 



AT if A/z = 0. 



(20) 
(21) 



For example, along a fixed boundary arc required by boundary 
conditions to be a line of constant h (e.g., entrance face OA 
of the dam) , nodes are adjusted so as to assume desired stream 
function values by using Eq. (21). 

The justification for Eqs. (18)- (21) is perhaps clear for situations 
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having no free surface, hence entailing no change in the (true) distri- 
butions h = h(x,y) and V = V(x,y) from iteration to iteration. In such 
cases the procedure shown here is nothing more than a Newton-Raphson 
iteration with a relaxation factor. A primary reason for the relaxation 
factor is, of course, to make it possible to prevent the free surface 
from shifting too abruptly and thus to insure that the current deriva- 
tives of h and V provide a meaningful basis for extrapolation. 

As a final point before turning to the subject of computer input, 
it should be noted that to avoid instabilities it may be necessary to 
use modified expressions for Ax and Ay at nodes where the derivative 
combinations appearing as denominators in Eqs . (19)- (21) tend to be 
small. This is the case, for example, along exposed discharge face BC 
of the dam when nodes are adjusted so as to assume prescribed stream 
function values. Here it is necessary to use 

or some similar expression, rather than Eq. (21), which ordinarily would 
be used in such a situation. (In theory, the denominator of Eq . (21) is 
identically zero at exit point B.) 

To describe actual computer input corresponding to the above pro- 
cedure, it is convenient to consider the special (but practical) case of 
a dam with no tailwater, i.e., a configuration in which downstream water 
level h 2 is reduced to zero, leaving the entire discharge face exposed 
to the atmosphere. A typical solution for such a configuration is illus- 
trated by Figure 12. Thirteen iterations were required to produce the 
final lattice of Figure 12c, at which point variations in the numerical 
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A 



5 ' z ' exit "^ max-' 




(a) Initial lattice 





(b) Lattice after 
first iteration 



(c] Lattice after 
final iteration 



Fig. 12. Graphical displays' illustrating the iterative gener- 
ation of a flow net for a darn with no tailwater . 



solutions for h and ¥ along the grid lines of constant i and j, respec- 
tively, were less than 0.2 percent. The appearance of the final lattice 
differs from that of a conventional flow net principally because the 
nodes are connected by straight lines rather than smooth curves and, in 
addition, because the grid lines do not represent equal increments of h 
and f. Here, in fact, the grid lines of constant j have ^-values which 
are distributed according to 



¥ = 



J 



.'-'max. 



Y/ 



t T 



Produced by "eps-p saved" on the ESL display console by typing 
'plot nil", then "plot 0$" (labels are added). See Section 7.2. 
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The lines of constant £ represent equal increments of h only in the re- 
gion between entrance face OA (i = 0) and the line of constant i whose 
upper end coincides with exit point B {i - i ex ±i) • For i ^ ^exit> ^~ 
values are determined by the elevations at which the lines of constant 
j intersect discharge face BC (see boundary conditions, Eq. (14)). 

To obtain a solution such as that of Figure 12, it is convenient 
to use the following procedure: 

(1) Specify lattice structure in i ,j -space (see Section 3.4), then 
define a mapping algorithm in terms of matrices "xx" and "yy" 
by typing 

def ine x=xx(i ,j) , y=yy(i,j) 

Note that it is easy to have the computer generate initial 
components for "xx" and "yy" (by means of appropriate do- 
loops) if, as in Figure 12a, the initial lattice is made up 
only of straight lines. 

(2) Specify governing equations and boundary conditions corre- 
sponding to Eqs. (13) -(16) using the boundary conditions 
h - y and Sf/Bn = 8j//94 along the estimated free surface. 
(Note that these conditions also apply along discharge face 
BC .) For example, with u-^ ■*■ h and 1*2 -+■ ¥, type 

set n=2, a(l,l)=a(2,2)=e(l,l)=e(2,2)=l, 
a(l,2)=a(2,l)=b=c=d=e(l,2)=e(2,l)= 
f=g»h=0 
impose 1 along 0,0, 0,1, Oj^ 
impose 2 along 0,^, 1,3 wax* ^max>° 
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do( do( do( set p (k, 1 ,m)=q (k, 1 ,m)=0 

for k=l step 1 unti 1 2) 

for 1=1 step 1 unti 1 2) 

for m=0 step 1 unti 1 2) 
set p(l,l,0)=q(l,l,l)=q(l,l,2)=l, 

q(2,2 ) 0)=p(2,2,l)=p(2,2,2)=l, 

r(l,0) = r(2,0) = r(2,l)=0, rd.l)^ 
define r (l ,2)=fy=f i t (y) , r(2,2)=ys 

where, as before, h-y is the upstream water level. (Note the 
use of the built-in function "ys", representing 3j//34.) 

(3) Obtain initial solutions for In. and ¥ by using form- and relax- 
commands as indicated in Section 6. 

(4) Define parameters "ph" and "pp" giving desired values of h and 
¥, respectively, in terms of the computed value "u (2 ,0 ,J max ) " 
of ^ at entrance point A, two interpolation functions "fi" and 
"fj", and the anticipated new elevations "yy (i ,£ max -i )+mu*dyy ( 
i,i' -i)" of nodes on discharge face BC by typing 

define ph=phf(i geq i exit ) , PP=fJ " u ( 2 >° >J m ax^ ' 

phf (0) =h ± +f i * (yy U exit , J max ) +mu"dyy U exit ,J max ) 
-hy) , phf (l)=yy ( i ,i max - i )+mu*dyy ( i ,i max - i ) 

Then specify "fi" and "fj" to indicate the desired contour 
distributions, e.g., type 

define fi = i/i exit . f J = 0/j max ) P° wer - 8 
to obtain the distribution used for the solution of Figure 12. 

(5) Specify parameters "dh" and "dp" representing the differences 
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between desired and computed values of h and V, respectively, 
by typing 

define dh=ph-u (l, i , j) , dp=pp-u(2, i , j) 

Then specify full corrections "dx" and "dy", as given by Eq. 
(19), for the interior and free-surface nodes by typing 

define dx=(hx*dh-hy*dp)/(hx*hx+hy*hy) , 
dy=(hy*dh+hx*dp)/(hx*hx+hy*hy) , 
hx=ux(l,i ,j) , hy=uy(l,i,j) 

Similarly, specify full corrections (i) "dxoa" and "dyoa", 
based on Eq. (21), for nodes on entrance face Ok, (ii) "dxbc" 
and "dybc", based on Eq. (22), for nodes on discharge face BC , 
and (iii) "dxoc", based on Eq. (20), for nodes on base OC by 
typing 

set CQ=xx(O,i max )/yy(0,j max ), 

cc= (xx U exit , J max ) -xx U max , 0) ) /yy (i exit , J max ) 

define dxoa=co*dyoa , dxoc=db/hx, 

dyoa=dp/(co*ux(2,0,j)+uy(2,0,j)) , 
dxbc=cc*dybc, dybc=dp*yy (i ,j)/u(2, i ,j) 

(6) Specify a macro "save" for computing all non-zero components 
of two correction matrices "dxx" and "dyy" by typing 

let save mean 
begin 
set i=0 
do( set dxx(0,j )=dxoa , dyy (0, j )=dyoa 

for j=l step 1 until j max -l) 
do( set i=£ max -j, dxx(i ,j)=dxbc, dyy ( i , j )=dybc 
for j=l step 1 until J max ) 
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set j=0 

do( set dxx( i ,0)=dxoc 

for i=l step 1 until i max - l) 

do( do( set dxx(i,j)=dx, dyy(i,j)=dy 

for i=l step 1 until i max ~J~l) 

for j=l step 1 unti 1 j max ) 
end 

Then specify the components of "dxx" and "dyy" which are to 
be zero for all iterations by typing 

set dxx(0,0)=dxx(7' max ,0) = 

dxx ( 'Jmax) =d Yy (0 'Jmax) ==0 
do( set dyy(i,0)=0 for i=0 step 1 until i max ) 

(7) Specify functions corresponding to Eq. (18) for calculating 
revised x ^-coordinates by typing 

def i ne newx=xx ( i , j )+mu-' c dxx( i , j ) 
newy=yy (i , j )+mu-'dyy (i , j ) 

Then, with a set -command, assign relaxation factor "mu" a 
starting value. Generally this value is not too critical as 
long as it is well below unity; 0.3 was used for the first 
few iterations of the solution shown in Figure 12. 

(8) Revise lattice geometry by typing 

save 

do( do( set xx (i ,j)=newx, yy(i,j)=newy 

for i=0 step 1 until i max ~}) 

for j=0 step 1 until i max ) 

and obtain new solutions for h and ^ by typing form- and 
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re lax- commands. 

(9) Repeat step (8), increasing "mu" occasionally, until adequate 
convergence is obtained. 

The following command specifies a recursive "adjust -macro" completely 
analogous to the one given previously with the procedure based on the 
use of interpolation formulae: 

let adjust (tol) mean 
begi n 

set i-£ exlt , j=«7 max st °P 
if abs(dybc) geq tol insert 

begin 

save 

do( do( set xx (i , j)=newx, yy(i,j)=newy 
for i=0 step 1 until "Z^^x-j) 
for j=0 step 1 until J max ) 

form relax adjust(tol) 

end 
end 

As with the previous procedure, convergence is generally slowest in the 
vicinity of exit point B. However, while errors are comparable in both 
procedures for a given number of nodes and for dams with steep discharge 
faces, they do not increase as rapidly with the present procedure when 
the slope of the discharge face is decreased. Further experiments are 
needed before more precise comparisons can be made. 

10.4. Over- and underrelaxation of iterated parameters. In the 
procedure for generating flow nets given above, a relaxation factor u 
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is needed to regulate the amounts by which the x,y -coordinates of the 
various nodes change from iteration to iteration. Such a factor may 
also be beneficial in other procedures, either to speed up convergence 
or simply to make convergence possible in cases where instability would 
otherwise be a problem. 

In general, a procedure which does not involve a relaxation factor 
can be modified so as to include such a factor as follows. Suppose, 
first, that the procedure normally uses a value (or set of values) a^ 
for some iterated quantity (or set of quantities) a in order to produce 
the -'I th solution. Then to introduce a relaxation factor y, it is neces- 
sary only to substitute for % per se the weighted average 

a^ = (l-y)a M + ya^ = a^ + y (a* - a^_ x ) . 

The modified procedure is identical to the normal procedure when y is 
one, generating in this case an increment Aa^ = (% - a^_]_) = (% - a^) 
in preparation for the -t th solution. However, for any single iteration 
the modified procedure will generate a smaller increment for a than the 
normal procedure if y is less than unity (underrelaxation) . Conversely, 
the modified procedure will generate a larger increment than the normal 
procedure if y is greater than unity (overrelaxation) . Depending on the 
situation, either of these effects may be desirable. For example, if a 
converges monotonically with the normal procedure, it may be possible to 
accelerate convergence by using overrelaxation. On the other hand, if a 
normally behaves in an erratic manner and/or diverges, underrelaxation 
may yield better results. In any event, the actual mechanics of intro- 
ducing a relaxation factor are quite straightforward. This can be illus 
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trated by reconsidering two of the examples presented earlier. 

A. Compressible flow. In the procedure for solving compressible 
flow problems (Example B, Section 10.2), a relaxation factor can easily 
be introduced at the point where the matrix components "vvm( i , j ) ", i.e., 
the current values of | grad $ | 2 at the various nodes, are recomputed. 
Specifically, it is necessary only to redefine the auxiliary parameter 
"vvfcn", used by macro "revise", in the following manner: 

define vvfcn=(l-mu)*vvm( i ,j )+mu- (ux(l , i ,j ) power 2 
+uy (l, i ,j)power 2) 

In test cases it has been observed that without a relaxation factor the 
"vvm(i,j) M tend to converge monotonically in much the same manner as in 
successive solutions obtained analytically by means of the Rayleigh- 
Janzen method. Therefore one would expect the optimum value of "mu" to 
be greater than unity. 

B. Nonlinear heat conduction. In the procedure for solving con- 
duction problems with temperature-dependent conductivities (Example A, 
Section 10.2), equation coefficients a^ = e-^_ = (u-|_/8 ) k Q are given 
directly in terms of the current solution matrix "u", which is overwrit- 
ten to produce the "next" solution during each iteration. Since using a 
relaxation factor to calculate some iterated quantity a^ requires the 
previous value a^_ 1; it is necessary to modify the original procedure so 
that previous values for either the coefficients a-Q = e-Q or the solu- 
tion u-i are saved. Choosing the former approach, the user could type 
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define a=e=f i t (k ( i , j ) ) , kfcn= (l-mu) *k ( i , j ) + 
mu*fc *(u(l,i ,j)/9 ) power s 

then have matrix "k" refreshed before requesting each new solution by 
using auxiliary parameter "kfcn" and a macro analogous to the revise- 
macro cited above. Clearly the added processing implied by this modi- 
fication must be weighed against the possible benefits of over- or under- 
relaxat ion. 

As a final point, note that it is not necessary to use the same 
relaxation factor with all iterated parameters. In fact, in cases where 
iterations serve more than one purpose, the added control provided by 
multiple factors may be very desirable. While no experiments have been 
performed to confirm this point, it seems likely that the procedure for 
generating flow nets which was discussed in Section 10.3 could be im- 
proved by using two factors, one to regulate the rate of change of the 
free-surface nodes and a second to control the changes of all other 
nodes . 



■112- 



11. Disk input and output 

EPS was designed strictly for on-line operation, but does provide 
a rudimentary mechanism for user-present disk I/O. With this mechanism 
it is possible to direct the computer (from a console) to take subse- 
quent input from a disk file prepared previously by standard CTSS proce- 
dures. Similarly, it is possible to cause the computer to produce a 
disk file which will, in essence, be a transcript of all subsequent con- 
sole interactions. Disk input is especially convenient in situations 
where several problems requiring common definitions are being treated, 
while disk output is useful when results generated by EPS are to be used 
as input to other programs. (Output files are written in a format which 
is compatible with standard CTSS editing and printing facilities, and 
generally can be made compatible with other application programs without 
difficulty. ) 

11.1. Disk input. Disk files which are to be used as EPS input 

must be written in line-marked format and must represent normal EPS com- 

f 
mands and macro calls with CTSS 6-bit character code. Such files can 

be created conveniently by means of the CTSS editing program, EDL, which 

is described thoroughly in [2], Sections AH. 3. 07 and AH. 9. 01. Two minor 

restrictions not specifically imposed by CTSS and EDL must be observed: 

(1) files must contain no lines with more than 83 characters, and (2) 

file names must not be reserved EPS vocabulary words (see Section 2) . 

Once a suitable input file, say name-^ name2, has been created, the user 



f 
See [2], sees. AB.2 and AC. 2. 
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can "resume eps" and, when it is appropriate, type 
read name-^ namp^^ 

to cause the computer to begin reading the file. Note that file names 
which normally would not be interpreted as identifiers by EPS must be 
introduced as literals. For example, to read the file "-cases 5a, b", 
type "read bed /--cases/ bed /5a,b/$". See Section 15.1. 

As with keyboard input, file reading proceeds until a line contain- 
ing a "$" is encountered, whereupon the input preceeding the "$" is pro- 
cessed. Normally this cycle is repeated until either all of the file is 
read or an error is detected, at which time control is returned to the 
user. An example is provided by the following sequence: 



edl 

W 936.7 
INPUT. 

let tabulated for v=b,d,e) mean 
beg i n 
i ndex 

expand begin v,z: end 
do(print v,z using format 

for v=b step d until e) 
i ndex 
end $ 
let index mean begin expand begin end end 
let format mean begin bed +f 4 . 2 ,8f9. 5*+ end$ 
define s i nh= (exp (x) -exp (-x) ) /2 , 
cosh=(exp (x)+exp(-x) )/2 
tabulate (s i nh , cosh for x=0,0.2,l)$ 

EDIT. 

file al pha beta 

R A. 050+2. 037 



] 



Cftzcutlon o{) 
input {iZz. 



resume eps 

W 945.2 
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PROCEED: 

read alpha beta$ 

FILE ALPHA BETA HAS BEEN OPENED. 

INDICATED SUBSTITUTION FOR 'TABULATE' WILL BE MADE. 

INDICATED SUBSTITUTION FOR 'INDEX' WILL BE MADE. 

INDICATED SUBSTITUTION FOR 'FORMAT' WILL BE MADE. 



X.SINH.COSH: 




0.00 


0.00000 


1.00000 


0.20 


0.20134 


1.02007 


0.40 


0.41075 


1.08107 


0.60 


0.63665 


1.18547 


0.80 


0.88811 


1.33743 


1.00 


1.17520 


1.54308 



Rtepon&e. 
to ■Input 
{ill. 



END OF FILE ENCOUNTERED. ALPHA BETA HAS BEEN CLOSED. 



PROCEED WITH CONSOLE INPUT: 
define tanh=sinh/cosh 
tabulate(tanh for x=.2,.2,l)$ 

X.TANH: 

0.20 0.19738 

0.40 0.37995 

ETC. 



Continu- 
ation oh 
nonmxt 
I/O. 



From this example it will be seen that EPS responds to commands in a 
disk file in exactly the same manner as it does to commands supplied 
directly from the keyboard. Thus, for most purposes, disk files can 
be considered logically equivalent to (parameterless) macros. In fact, 
they offer two advantages over macros: (1) they do not consume core 
storage and (2), since they can be prepared with a text editor, they 
are easier to specify. Of course, they also have a disadvantage, name- 
ly, increased processing time, and generally should not be used to 
represent input sequences that are to be repeated often. 

One further point should be mentioned. When an error is detected 
while input is being taken from a disk file, EPS prints whatever error 
message it would normally print, then returns to the user's console 
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for further instructions. It does not, however, close the input file 
in such a situation, but, in fact, will continue reading the file from 
where it left off if given the command 

read swi tch$ 

Thus it is frequently possible to correct the source of an error, then 
proceed with the reading of an input file, without serious loss. 

11.2. Disk output. EPS does not provide disk output as an alter- 
native to console output, but simply offers a special I/O mode in which 
all information printed at the user's console, including input typed by 
the user himself, is written concurrently into a disk file. To enter 
this special mode, the user types a command of the form 



record name name $ 

where name name identifies the file to be written. Provided the user's 
track quota is not exhausted, the computer will respond by printing 



Actually, the command "read switch$" has more general utility. 

Disk files are allowed to be either active, i.e., in use, or inactive, 

and the command "read switch$" simply indicates that the present status 

of an input file is to be reversed. Thus it is possible not only to 

reactivate an input file which has been deactivated because of an error, 

but also, by placing a "read switch$" in an input file, to have the 

file deactivate itself in order to allow for user interaction, e.g., 

specification of parameters needed for further processing. Of course, 

after such an interaction, the user must issue a "read switch$" command 

himself in order to return control to the disk file. 

§ 
If no name* name„ exists, one will be started. Otherwise, 

output will be appended to the existing name* name,*. Hence, it is 

possible to accumulate output from several sessions in the same file. 

See [2], sec. AH. 1.04. 

-116- 



FILE NAME NAME HAS BEEN OPENED. 

After this, all information typed by the user's console printer will 
be copied onto the disk, except as noted below. 

(1) Input cancelled by means of the CTSS "erase" characters 
(" and #) and "kill" characters (? and @) never reaches 
EPS and consequently is not recorded (see Section 2). 

(2) The recording of all I/O can be discontinued by typing 

record switch$ 

and restarted by typing the same command again. Thus the 
recording of interactions which are not of interest can 
be suppressed at will (see footnote on preceeding page). 

Files created by use of the record-command are line-marked and written 

in CTSS 6-bit character code. Therefore, they can be printed by means 

t 
of the CTSS PRINT- command and edited with EDL. 



+ See [2], sec. AH. 5.03. 
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12. i n[»ut i'tr i fi cat i on 

In ir.ost applications it becomes desirable from time to time to 
have the computer report its current status with respect to a particu- 
lar class of input or with respect to a given item within a particular 
class. In extended or interrupted sessions, for example, it may be 
helpful simply to be reminded of what macro identifiers have been intro- 
duced in past input. Or, as a precaution against typing and transmis- 
sion errors, it may be useful to have the computer report back compli- 
cated curve descriptions or parameter definitions. Operations such as 
these can be carried out by using the EPS list- and review-commands, 
as discussed below. 

12.1. Determining active identifiers. It is possible to obtain 
listings of all active identifiers in each of one or more classes c-^, 
<?2> ■•-, <?//, where the possible classes are "curves", "variables", 
"macros" and "commands", by typing 



i st cj 1 , c 2 , ..., c N 



For each c»f, for which active identifiers exist, the computer will print 
a message of the form 

C\ CURRENTLY ACTIVE: 'NAME\, NAME' 2 , ■■■, NAME?/. 

Thus, for example, to determine what curve and variable identifiers are 
currently in use, type 

list curves, variables $ 
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A possible response would be 



CURVES CURRENTLY ACTIVE: 'ALPHA, BETA 1 . 

VARIABLES CURRENTLY ACTIVE: 'LIMIT, DELTA, OMEGA, 
R, Q, P, H, G, F, E, D, C, B, A, Y, X, J, I, U, UX , 
UY, FLUX, UXX, UYY, UYX , UXY, FIT, XS, YS , OCTANT'. 



(The word "variables" is something of a misnomer, for its use causes 
constant parameters and built-in functions, as well as variable param- 
eters, to be listed.) Identifiers that have been introduced by the user 
are always listed in reverse chronological order, the first identifier 
listed being the one most recently introduced, and so forth. 

12.2. Determining individual definitions. It is possible to ob- 
tain output giving current definitions for one or more curves, variables, 
or macros by typing a command of the form 



rev i ew 



curve 


ar i able 


macro 



begin dj_, dj, . . . , d^ end 



(The braces Q are not meant literally, but merely signify that a choice 
is available.) The dfc may be single identifiers or they may be identi- 
fiers followed by certain qualifiers, depending on what class of infor- 
mation is being reviewed and on what output is desired. In particular, 
when curves are being reviewed, the dfo may assume two forms: 



cuvvename or cuvvename along i]_,j]_, ^2"? 2' ^3»<7 3 

The second, or qualified, form is used when output is to be restricted 
to the portion of cuvvename passing from i\,3\ through i 2 j J 2 t0 ^3>^3- 
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Similarly, when variables are being reviewed, the dfc may take the form 
pararrmame or pararrmame (s^ ,s 2 , • • • ,s M ) 

where the s^ represent subscripts and may, in general, be arbitrary alge- 
braic expressions. It should be noted that the number of subscripts M 
need not be as great as the actual number of subscripts associated with 
a given parameter. In fact, by omitting one or more subscripts, the 
user can obtain output indicating the current dimensions of a subscripted 
quantity (see example below). Finally, when macros are being reviewed, 
each d k must be a single identifier; argument lists, for example, are not 
permitted here as they are with the expand -command discussed in Section 

9.5. 

Output generated by means of the review-command is largely self- 
explanatory, as is illustrated by the following sequence: 

PROCEED: 

let mtx(w=a,b,c,d ,) mean 

beg i n 

set w(l,l)=a, w(l,2)=b, w(2,l)=c, w(2,2)=d 

end$ 
INDICATED SUBSTITUTION FOR 'MTX' WILL BE MADE. 

PROCEED: 

review macro begin mtx end$ 

MTX(W=A,B,C,D,) MEANS SET W(l,l)=A, W(l,2)=B, W(2,l)=C, 
W(2,2)=D 

PROCEED: 

set n=2, b=c=d=f=h(2)=0 

mtx (a=l, 0,0,-1, e=l, 0,0,-1, g=0, 0,1,0) 

define fx=fit(x), h(l)=po*(3.75+fx*fx)/s 

review variable begin a ,b ,e ,g ,h (l) ,h (2) end$ 

A=(A(1,1)...A(2,2)) 

B=0 

E=(E(1,1)...E(2,2)) 
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G-(G(1,1)...G(2,2)) 

H(l)=PO*(3.75+FX*FX)/S 

H(2)=0 

PROCEED: 

review variable begin a (1,1) ,a(2,2) ,e(l) ,e(2) ,po,fx,s end$ 

A(l,l)-1 

A(2,2)— 1 

E(1)-(E(1,1)...E(1,2)) 

E(2)=(E(2,1)...E(2,2)) 

■PO' NOT SPECIFIED. 

FX=FIT(X) 

'S' NOT SPECIFIED. 

PROCEED: 

append 0,0, 3,0, 3,3, 0,3 to ca close ca 

append 1,1, 2,1, 2,2, 1,2 to cb close cb 

impose k along 0,0, 2,0, 3,0 

impose 5 along 0,0, 0,2, 3,0 

impose 6 along 1,1, 2,2, 1,1 

impose 7 at 3,0 $ 

DEFINITION OF NEW CURVE 'CA' HAS BEEN COMPLETED. 

'CA' HAS BEEN CLOSED. 

DEFINITION OF NEW CURVE 'CB 1 HAS BEEN COMPLETED. 

'CB' HAS BEEN CLOSED. 

PROCEED: 

review curve begin ca along 0,1, 0,0, 3,2, cb end$ 

GRID COORDINATES AND LOADING INDICES OF 'CA': 
( 0, 1) ( 5, 5) 



( o, 


0) 


( 


5, k) 


( 1, 


0) 


( 


A, k) 


( 2, 


0) 


( 


<t, h) 


( 3, 


0) 


( 


7, 7) 


( 3, 


1) 


( 


5, 5) 


( 3, 


2) 


( 


5, 5) 


RID COORDINATES AND LOADING INDICES 


( 1, 


1) 


( 


6, 6) 


( 2, 


1) 


( 


6, 6) 


( 2, 


2) 


( 


6, 6) 


( 1, 


2) 


( 


6, 6) 



CB' 



It will be seen that curve reviews give two segment numbers (i.e., 
"LOADING INDICES") for each boundary node. The first segment number 
corresponds to the boundary arc immediately preceeding the node, and 



-121- 



the second to the arc immediately following the node, where direction of 
travel is implied by the order in which the nodes are listed. The curves 
specified in this example would not be specified in practice, since they 
correspond to a lattice in which all nodes are boundary nodes! However, 
the equation parameters have practical significance: they correspond to 
the plate equation with constant bending stiffness "s" and variable pres- 
sure "po- c (3 • 75+fx ; -fx) " (see Section 4.3, Example C) . 
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13. Program interruption and run termination 

Occasionally the user will wish to prevent or curtail the process- 
ing of a unit of input, either because he recognizes that he has made an 
error or simply because he changes his mind. Further, he will wish at 
some point to quit EPS completely and have control passed back to the 
CTSS supervisor. These operations can be accomplished by means of the 
program interrupt button and the EPS quit-command, as described below. 

13 . 1 . Interrupts . Normal processing can be discontinued at 
any time by pressing once the program interrupt button (labeled "INT", 
"BREAK", "RESET", or "INTERRUPT", depending on console type). This 
action causes the computer to ignore any commands supplied since the 
last "PROCEED:" but not yet executed, and with certain minor exceptions 
results in the response 

INT. 
PROCEED: 

After this, normal operation of the program can be continued. It is 
emphasized that interrupts have no ability to "undo" commands once they 
have been executed. Thus, in particular, data specified in previous 
interactions is neither lost nor disturbed in any way. 

The user should note that if the interrupt is pressed during 
the specification of input, i.e., before a "$" and final carriage re- 
turn have been given, no commands supplied since the last "PROCEED:" 
will have been executed. Thus in this situation the interrupt is anal- 
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ogous to the CTSS "kill" characters, ? and @ (see Section 2), except 
that it causes all preceeding lines of the unit of input, as well as 
the current line, to be ignored. (Also, of course, it causes output 
to be produced, including a new "PROCEED:". The latter must be awaited 
before typing can be resumed.) 

On the other hand, if the interrupt is pressed after a "$" and 
final carriage return have been given, any number of commands, includ- 
ing none, may in general have been executed. In most cases, output 
produced before the interrupt will indicate how far processing has 
gone; if not, appropriate print- and review-commands may be used to 
ascertain this information (see Sections 7 and 12). 

To avoid confusion, the user should be warned that interrupts 
often do not "take effect" immediately. Typically, the machine will 
require several seconds, say 5 to 15, to respond with "INT. 0", then 
possibly several more to produce a "PROCEED:". Moreover, because EPS 
normally completes the command it is currently executing (if any) before 
actually obeying an interrupt, it is possible to receive a small amount 

J 

of regular output after the line "INT. 0" and before the new "PROCEED:". 
Therefore, patience is necessary; attempting to hurry the machine by 
pressing the interrupt repeatedly will either delay matters further 
or, as indicated in the next section, produce a completely different 
effect . 



J EPS allows itself to be interrupted in "mid-command" only when 
it is generating or solving finite-difference equations. See descrip- 
tions of the form- and re 1 ax-commands in sec. 15 for further details. 
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13.2. Run termination. Control can be returned to the CTSS 
supervisor at any time by pressing the interrupt button twice , or, at 
times when EPS is in input status, by typing the command 

quit $ 

It is necessary to use the interrupt button if immediate termination 
is desired at a time when an extended calculation is in progress. How- 
ever, for reasons given later, it is recommended that the quit-command 
be used whenever possible. In either case, the computer will ultimately 
respond with a line of the form 

R AT-1+AT2 

where A^ and NF 2 8 ive computation time and overhead, respectively, in 
seconds. After this, commands must be directed to CTSS. For complete- 
ness, three common situations are considered. 

(1) The user may realize immediately after quiting EPS that he 
has not finished his computations. To continue, he can type 

start 

whereupon the computer will resume execution of EPS at the 
point where it left off. Note that this procedure can be 
used deliberately to determine how much computer time is be- 
ing used. 

(2) The user may anticipate continuing his computations at a 



+ Cf. [2], sec. AH. 7. 03. 
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later time, in which case he can type 



save name-^ 

to create a disk file "name^ saved" to be resumed in the 
future. In applications requiring several iterative 
stages, it may be wise to create "saved" files at various 
points along the way in order to avoid having to recalcu- 
late earlier results if later strategies fail. 

(3) The user may be satisfied with his results, certain that 
he is finished, and ready to leave his console, in which 
case he can log out as per [2], Section AH. 1.02. 

It should be noted that output files opened by means of the record- 
command (Section 11.2) are closed when and only when a quit-command 
is received; such files may be lost if the interrupt button is used 
to terminate a run. Further, if EPS is restarted (or saved and re- 
sumed) after being terminated by means of the interrupt button, it 
may enter input status without printing the usual "PROCEED:". To 
determine whether or not the computer is, in fact, awaiting input 
in this and other situations, type a "$" and carriage return. This 
has no computation effect, but insures that the computer will print 
a "PROCEED:" when it next needs input. 



+ Cf. [2] , sec. AH. 3. 03. 



Disk files created in this manner may be from approximately 
30 to 78 records in length, depending on which version of EPS is be- 
ing used and on how much information has accumulated in core memory 
as a result of previous interactions. 
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14. lirror messages 

Whenever continued processing of an input sequence is impossible 
because of bad form, invalid or insufficient data, or the detection of 
an illegal machine condition, the computer prints an error message indi- 
cating specifically what difficulty has been encountered. Depending on 
what stage processing has reached and on the particular difficulty that 
has arisen, error messages may either demand a special response from the 
user or simply state that processing has been cut short and ask, by way 
of a "PROCEED:" message, for more input. In either case, messages are 
in principle completely self-explanatory, except in cases of input con- 
version errors, which are reported in terms of a numerical error code // 
and a line image L as follows: 



ERROR N DETECTED IN SCANNING: 

L 

RETYPE THIS LINE AND ANY OTHERS WHICH FOLLOW: 



Usually messages of this type are incurred as a result of poor console 
technique ( see Section 2), as may be suggested by the list of possible 
codes given below. 



ERROR CONDITION 



Memory overflow (very unlikely in this context) 

Incomplete or possibly truncated word or number 
(usually a line with more than 83 characters) 

Number greater than 1.701412E+38 

Illegal character sequence, e.g., "3^beta" 



Such messages are rare once the user has gained a little experience. 
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Much more likely, in general, are error messages resulting from 
the inadvertent omission of elements in algebraic expressions. Fortu- 
nately, the computer is more generous in such cases. All algebraic in- 
put to EPS is checked for syntax errors, i.e., missing (or excessive) 
parentheses, operators, and operands, and "corrected" when necessary 
by means of a simple but surprisingly effective strategy. "Corrected" 
expressions are reported with messages of the form 

INCORRECT INPUT STRING MODIFIED TO READ 'EXPR'. 
IS THIS WHAT YOU MEAN. . . 

in response to which the user is expected to type "yes" or "no". If 
he types "yes", the computer will resume its normal operation, so pre- 
sumably no loss will be suffered. If he types "no", the computer will 
execute a fatal error sequence , which means it will discard the command 
containing the erroneous expression, as well as any subsequent commands 
waiting to be processed. This action may or may not be explicitly ac- 
knowledged in further error messages, but in all cases it will neces- 
sarily take place. Note that it is perfectly legitimate to answer 
"yes", even when an improper correction has been made, in order to avoid 
retyping valid input which would otherwise be discarded. Afterwards 
improperly corrected information can be respecified, as shown in Sec- 
tion 8. 

Fatal error sequences are executed in general when commands are 
found to contain non-correctable syntax errors, and also when they are 
given in terms of invalid data. However, non-fatal error sequences are 
normally executed when commands cannot be completed because of insuf- 
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ficient data, inadequate memory, or any other conflict potentially re- 
solvable through subsequent user interaction. In such situations, input 
sequences of more than one element which have not been completely proc- 
essed are not discarded, but instead are saved in a special push-down 
list called the command stack. Processing of the latest, i.e., top- 
most, input sequence on the command stack can be reinitiated at any time 
by means of the command "retry", while deletion of the top-most sequence 
can be effected by means of the command "ignore". In practice this ar- 
rangement is both natural and convenient, as the following example may 
illustrate: 



PROCEED: 

define radi us=i*sqrt (c2+sqrt (c2*c2+15) ) /10 

c2=cos(2*t), t=j*pi/20 
set i=10, j=0 print t, radius 
set j = 10 print t,radius$ 

PROCESSING HALTED ON CALL FOR UNDEFINED QUANTITY 
'PI 1 , WHICH WAS REFERENCED BY 'T 1 . 

REDUCTION OF 'T, RADIUS' CANNOT BE COMPLETED. 

COMMAND STRING 'PRINT T, RADIUS SET... 1 HAS BEEN 
SAVED. TYPE 'RETRY' TO REINITIATE PROCESSING, 
'IGNORE' TO HAVE STRING DELETED. 



PROCEED: 

set pi=A*atan (1) retry$ 

0.00000000000 2.2360683E+00 

1.5707963E+00 1.7320509E+00 



Clearly, command sequences which cause non-fatal errors will accumu- 
late indefinitely in the command stack unless they are successfully 
retried or deleted by means of an "ignore". To avoid unnecessary sto- 
rage demands, the latter command should be typed immediately if it is 
recognized that a given sequence will not be retried. 
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15. Command descri ptions 

In this section, descriptions of all of the EPS commands are pre- 
sented in outline form. These descriptions take a broader view of the 
system than does the material of preceeding sections. In addition, they 
introduce a certain amount of new information, e.g., optional input 
forms, brief explanations of how the commands are actually processed, 
and listings of the conditions under which error sequences are executed. 
Thus this section should be of particular interest to experienced users 
who wish to pursue new applications. Furthermore, it brings together 
in a form convenient for quick reference information which may be useful 
to any user who encounters difficulties while using the computer. 

15.1. Terminology. For each command, descriptive material is en- 
tered in a fixed order under the following headings: 

PURPOSE NORMAL RESPONSE 

REFERENCES SIDE EFFECTS 

FORMAT INTERRUPTS 

EXAMPLES NON-FATAL ERRORS 

NORMAL OPERATION FATAL ERRORS 

Entries are omitted when inappropriate, e.g., if a command normally pro- 
duces no response, the NORMAL RESPONSE entry is omitted. 

The material is largely self-explanatory, except possibly for some 
of the notation and terminology used in certain FORMAT entries. For 
most commands, the FORMAT entries define admissible input forms in two 
steps. First, with a notation which is a slight extension of that em- 
ployed in preceeding sections, basic structure is defined by an expres- 
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si on in which: 

(1) information intended literally is shown in sans-serif type; 

(2) command variables, or arguments, are shown in italics or 
other easily distinguished type styles; 

(3) alternatives are listed vertically between braces Q ; and 

(4) information which may be omitted under certain conditions 
is enclosed in square brackets [] . 

Thus, for example, the FORMAT entry for the delete-command begins with 
the expression 



delete 



[aurvename] 
H»j'l> l>2»<72»] ^3»<?3 [ from aurvename] t 



This indicates that a delete-command may assume the following basic 
forms : 

delete 

delete aurvename 

delete ii,j\, 13,^3 

delete i\,j\, ^3»J3 from aurvename 

delete i-\_,Ji_, i2>3l> ^3»J*3 

delete i]_,«7]_, ^2»^2» ^3 » J3 from aurvename 

The appropriateness of any given form depends, of course, on the action 
desired by the user, as well as on conditions explained in the NORMAL 
OPERATION entry. 

After basic structure has been defined, admissible forms of argu- 
ments, where allowed, are specified verbally. Thus in the case of the 
delete-command, the FORMAT entry concludes with the phrase 
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where curvename is an identifier and the ifc,3\i are algebraic 
expressions (1 < fc < 3) . 

The meanings of terms such as "identifier" and "algebraic expression" 
may at this point be intuitively clear, but for completeness are given 
informally below, along with various related definitions. Starred items 
define terms actually used in command descriptions. 

letter: one of the 26 letters of the alphabet, a single quota- 
tion mark ('), a form feed, or a period. 

digit: one of the 10 decimal digits. 

printed delimiter: one of the following: +-/*()=: , $ 

character: a letter, digit, printed delimiter, blank, or tab. 

normal word: a string of up to 83 letters and digits beginning 
with a letter. (Restriction: if the string begins with a period, 
it must consist of only the period, or else the period must be 
followed by a letter.) 

literal: a string s of up to 77 characters introduced by means 
of an expression of the form "bed As A", where A is any character 
not present in s and where the space between "bed" and the first 
A must be a single blank, not a tab or carriage return. 

*word: a normal word or literal. 

^integer: a string of one or more digits. 

number with decimal point: an integer followed by a decimal point 
(i.e., period), a decimal point followed by an integer, or an in- 
teger followed by a decimal point followed by an integer. 

number with scale factor: an integer or number with decimal point 
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followed by an "e" followed by a sign (+ or -) followed by an 
integer. The sign may be omitted when "+" is intended. 

number: an integer, number with decimal point, or number with 
scale factor. 

* input element : a word, number, or printed delimiter. 

*input sequence: a sequence containing any number of input ele- 
ments, including none. (Blanks, tabs, and/or carriage returns 
must be used to separate adjacent words and numbers, but are 
not considered part of the sequence per se.) 

reserved command word: any word shown in sans-serif type in the 
summary of commands given in Table 1, Section 16. 

function name: any word shown in sans-serif type in the summary 
of built-in functions given in Table 2, Section 16. 

multiliteral operator : any word shown in sans-serif type in the 
summary of operators given in Table 3, Section 16. 

reserved word: a reserved command word, function name, or multi- 
literal operator; also any one of the following: if insert bed 

* identifier: any word which is not a reserved word. 

* parameter designator : an identifier, or, if the parameter is sub- 
scripted, an identifier followed by an expression list enclosed 
in parentheses. 

function call : a function name, or, if the function requires ar- 
guments, a function name followed by an expression list enclosed 
in parentheses [see Table 2, Section 16). 

unary operator : any input element designated as a unary operator 
in Table 3, Section 16. 
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diadic operator : any input element designated as a diadic opera- 
tor in Table 3, Section 16. 

* algebraic expression : a number, parameter designator, or function 
call, or, proceeding recursively, a sequence of the form 

o 1 e or e o 2 e ' or (e) 

where o l is a unary operator, o 2 is a diadic operator, and each 
of the e's is itself a legitimate algebraic expression. For oper- 
ator precedence information, see Table 3, Section 16. 

expression I ist : a single algebraic expression or a sequence of 
algebraic expressions, each separated from the next by a comma. 



15.2. Descriptions. Command descriptions follow. Underlined 
section numbers listed as REFERENCES indicate sources of background in- 
formation and/or examples which are particularly relevant. The headings 
FATAL ERRORS and NON-FATAL ERRORS refer specifically to conditions under 
which fatal sequences and non-fatal error sequences, respectively, are 
initiated (see Section 14) . 
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append -command 

PURPOSE: To define or extend a bounding arc. 
REFERENCES: Section 3, especially Z_A, and Section 8.1. 

FORMAT: append ii,j\, i2'3l> •••» ^N'JN t to cuvvename] 

where the "£fc,Jfc are algebraic expressions [1 - k - N, N ^ 1) and 
where cuvvename is an identifier. 

EXAMPLES: append "5,0, 5,0, 5,9, "5,9 to zeta 

append 0,0, imax,0, imax+jmax, jmax, 0,jmax 

NORMAL OPERATION: The ifcjjfc are evaluated, rounded to nearest integers, 
and interpreted as the grid coordinates of successive vertices (in i ,j - 
space) of a bounding arc. A list of all nodes on the arc, including 
nodes on the grid lines and grid diagonals formed by connecting the ver- 
tices, is constructed. If the phrase "to cuvvename" is specified but 
identifier cuvvename is not presently- active, the node list is taken as 
the definition of a new arc. If "to cuvvename" is specified and found 
to refer to an existing, unclosed arc with an end point which coincides 
with either i\,3\ or i^,jpj, then the existing arc is extended in accord- 
ance with the new node list. If "to cuvvename" is omitted, a new arc 
with a computer-generated name is defined. 

NORMAL RESPONSE: Depending on whether a new definition or an extension 
is introduced, either 

DEFINITION OF NEW CURVE 'CURVENAME' HAS BEEN COMPLETED, 
or 

ADDITION TO 'CURVENAME' HAS BEEN COMPLETED. 

SIDE EFFECTS: Deletion of the difference equation coefficients and 
solution matrix generated as a result of any previous form-command and 
of the data base generated as a result of any previous plot-command. 

NON-FATAL ERRORS: Memory overflow and references by the £fe,Jfe to unde- 



■135- 



append-command (continued) 



fined quantities. 



FATAL HRRORS: Non-correctable syntax errors (e.g., only one pair of 
-^/-coordinates) and illegal data (e.g., successive vertices which can- 
not be connected by a single grid line or grid diagonal). 
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close -command 

PURPOSE: To cause the end points of a bounding arc to be connected so 
as to form a closed contour. 

REFERENCES: Section 3, especially 3.4 , and Section 8.1. 

FORMAT: close [eurvename] 

where eurvename is an identifier. 

EXAMPLES: close zeta 
close 

NORMAL OPERATION: The end points of contour eurvename, if eurvename is 
specified, or of the only contour which exists (presuming only one ex- 
ists), if eurvename is omitted, are inspected. If these points are dif- 
ferent and can be joined by a grid line or grid diagonal, then an appro- 
priate node list is generated and inserted between the existing end 
points. If the existing end points are identical, they are simply coa- 
lesced. 

NORMAL RESPONSE: 'CURVENAME' HAS BEEN CLOSED. 

NON-FATAL ERRORS: Memory overflow and attempts to close contours whose 
end points cannot be connected by a single grid line or grid diagonal. 

FATAL ERRORS: Omission of the argument eurvename when more than one 
contour exists and specification of a eurvename which does not exist. 
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def i ne-command 

PURPOSE: To specify functional relationships between algebraic param- 
eters . 

REFERENCES: Sections 3^, k_/2_, 4.3 , 5.3 , 7.1, EL2, 10.2, and 10.3. 

FORMAT: define p 1± = p 12 = ■■■ = Pm ± = e l> ?21 = ?22 = ••• = 
P2M 2 = e 2 ' ■ ' ' ' PN1 = PN2 = • • ' = PNM N = e N 
where the pun are parameter designators and the e^ are algebraic 
expressions (1 < k < N, N > 1; 1 < I < i^, M^ > 1). 

EXAMPLE: define a=e=2+fx*fx, fx=fit(x), r (2 ,m)=fx/3 .^7 

NORMAL OPERATION: Argument groups of the form 

Pkl = Pkl = ••■ = Pkl4 k = e k ^ e -g- "a=e=2 + fx>vfx") 

are rewritten so as to become 

PkM k = e k> PkM k -l = PkM k > ■■•'Pkl- Pkl ^ e -g- "e=2+fx*fx, a=e") 

Thus, in effect, the original command is made to assume the form 

define p 1 = e ] _, p 2 = e 2> ■> PL = e L 

For each assignment p^ = e^, a subroutine for calculating e^ is compiled 
and the data base arranged so that in the future this subroutine is exe- 
cuted each time a value for p^ is required. If p^ has subscripts, these 
are evaluated and e^ is associated with the particular component that is 
indicated. Thus if "m" is 3 when the command of the example is speci- 
fied, "fx/3.47" is associated with the single matrix component "r(2,3) M - 

SIDE EFFECTS: Previous definitions are deleted and subscript ranges are 
modified as needed. 

NON-FATAL ERRORS: Memory overflow and references in subscripts of the 
p's to undefined quantities. 

FATAL ERRORS: Non-correctable syntax errors. 
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delete-command 

PURPOSE: To delete all or part of a bounding contour. 
REFERENCES: Sections 3 and 8.1. 



FORMAT: delete 



[curvename \ 
>l»Jl» l>2>^2'] ^3 "73 [ from curvename] 



where curvename is an identifier and the t\i,3\i are algebraic 
expressions (1 < fe < 3) . 

EXAMPLES: delete 

delete zeta 

delete 0,5, 0,6 from curve2 

delete 0,jmax, 2,jmax, imax,jmax 

NORMAL OPERATION: If no argument is specified and only one contour ex- 
ists, that entire contour is deleted. If the single argument curvename 
is specified, the entire contour curvename is deleted. If a phrase of 
the form "ii,j\, 'i-2i32> ^3»t?3 f rom curvename" is specified, the £fc,Jfc 
are evaluated and rounded to nearest integers, then the portion of 
curvename running from node ii,Ji through -i 2 » J 2 to "^ 3 > J 3 is deleted. 
The i-\i,d\i may specify an arbitrary arc if the contour in question is 
closed; however, if the contour is unclosed, the specified arc and the 
unclosed contour must have at least one end point in common--i .e. , it is 
not permissible to delete an intermediate arc from an unclosed contour, 
breaking the contour into two pieces. The phrase "i2>>72 j " can ^ e omit- 
ted at any time if the contour is unclosed, but can be omitted when the 
contour is closed only if end points i\,3\ and ^3,^3 are neighboring 
nodes (in which case the single line connecting these nodes is deleted) . 
The phrase "from curvename" can be omitted whenever only one contour 
containing the specified arc exists--i . e. , almost always, in practice. 

NORMAL RESPONSE: Depending on whether an entire contour or a portion of 
a contour is deleted, either 
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de 1 ete-command (continued) 



'C!JF7t:iJAME' HAS BEEN DELETED, 
or 



DELETION FROM 'CUR7FNAME' HAS BEEN COMPLETED. 

SIDE EFFECTS: Deletion of the difference equation coefficients and 
solution matrix generated as a result of any previous form- command, of 
the data base generated as a result of any previous plot-command, and 
of the segment numbers associated with all deleted arcs. 

NON-FATAL ERRORS: References among the ifc.Jfe to undefined quantities. 

FATAL ERRORS: Omission of arguments when more than one contour exists 
and specification of invalid or ambiguous data. 
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expand-command 

PURPOSE: To test macro definitions and cause printing of arbitrary text 
strings (e.g., table headings, etc.). 

REFERENCE: Section 9, especially 9.5 . 

FORMAT: expand begin p end 

where p is any input sequence which does not contain a "$" or un- 
matched begin-end pairs. 

EXAMPLE: expand begin this is a text string, end 

NORMAL OPERATION: The sequence p is passed through the EPS macro pro- 
cessor so that any macro calls it might contain are transformed into 
corresponding substitutions. 

NORMAL RESPONSE: The sequence resulting from the expansion of p, if 
any, is printed. If p is the null sequence or "expands to the null se- 
quence" (see comment -macro, Section 9.2), a blank line is printed. 

FATAL ERRORS: Memory overflow and references to undefined quantities in 
the determining expressions of if -phrases (see Section 9.3). 
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form-command 

PURPOSE: To generate finite-difference equations. 

REFERENCES: Section 6, especially 6.1 , and Section 8. 

FORMAT: form 

NORMAL OPERATION: Existing difference equation coefficients, if any, 
are discarded and a tally-command is simulated. If no solution matrix 
exists (necessarily the case if the form-command is being used for the 
first time or if an append- or delete-command has been given since the 
form-command was last used), one is created and set to zero. Similarly, 
if a solution matrix does exist, but the number of dependent variables 
"n" has been changed, the old matrix is discarded and a new one is cre- 
ated and set to zero. Otherwise, the existing matrix is left intact. 
Then (new) difference equation coefficients for the current lattice and 
differential system are generated. Coefficients corresponding to all 
"n" partial differential equations or boundary conditions are produced 
at first one node, then another. In terms of the lattice configuration 
in i, j -space, processing is from left to right along each grid line of 
constant j, starting with the line of smallest j, say j m ± n , and proceed- 
ing to ■7'min +1 > <7min+ 2 > etc - 

NORMAL RESPONSE: The response, if any, resulting from the simulated 
tally-command, then, if a (new) solution matrix is created, the line 

SPACE FOR SOLUTION MATRIX HAS BEEN ALLOTTED AND ZEROED. 

INTERRUPTS: The user can give a single interrupt signal (see Section 
13.1) either to curtail processing entirely, as when he discovers a mis- 
take, or simply to determine how far along processing has come. The 
computer will respond by printing 

GENERATION OF DIFFERENCE EQUATIONS HALTED AT GRID POINT 
(I, J). DO YOU WISH EQUATION GENERATION TO CONTINUE... 

after which the user must type "yes" or "no". 
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form- command (continued) 

NON-FATAL ERRORS: Illegal conditions recognized by the tally-command, 
memory overflows, and missing or unacceptable data (e.g., lattice param- 
eters "x" and "y" which describe a lattice with illegal geometry) . 
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ignore -command 

PURPOSE: To discard a command sequence saved as a result of a non-fatal 
error. 

REFERENCE: Section 14. 

FORMAT: ignore 

NORMAL OPERATION: I'hc most recently saved command sequence is discarded 
and the one proceeding that, if any, is readied for possible retrial. 

FATAL ERRORS: i g no re -commands given when no saved command sequences 
exist to be ignored. 
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i mpose-command 

PURPOSE: To associate a segment number, hence a given set of boundary 
condition parameters, with a bounding arc. 

REFERENCES: Section 5, especially 5.2 , and Section 8.3 . 



FORMAT: impose m 



[of curvename] 



along i^Ji, M-^^lA ^3»J3 
at ii,j\ 

where m and the ifc,jfc (1 ^ k ^ 3) are algebraic expressions and 
ourvename is an identifier. 

EXAMPLES: impose 3 along 0,0, 2,0, imax+jmax, jmax 
impose newm at 5,9 of zeta 

NORMAL OPERATION: The value of expression m is obtained, rounded to the 
nearest integer, and interpreted as a segment number to be associated 
with the arc specified by the remaining portion of the command. If the 
remaining portion has the form "along 1^,3 \, ^2"7*2' ^3»<73 °* euvvename" , 
the ifc,j \ are evaluated and rounded to nearest integers, and m is asso- 
ciated with the portion of ourvename running from i\,3\ through ^2 jt7 '2 t0 
t 3 ,j 3 . The phrase ""Z-2,J2>" can be omitted at any time if the contour in 
question is unclosed; however, if the contour is closed, it can be omit- 
ted only when i\,3\ and £3, J 3 are neighboring boundary nodes (in which 
case m is applied only along the line connecting the two neighbors) . If 
the remaining portion has the form "at i\,3\ of ourvename", i± and q\ 
are evaluated and rounded to nearest integers, and m is associated with 
the arc which includes ^i,3\ and extends half-way to each of its neigh- 
bors (see Section 5.3, remark (3)). In all cases, the phrase "of euvve- 
name" can be omitted whenever the arc being specified exists on only one 
bounding contour- -i .e. , the phrase can be omitted almost always in prac- 
tice. 

SIDE EFFECTS: Existing segment numbers are overwritten without comment. 

NON- FATAL ERRORS: References among the ifc,j\ to undefined quantities. 
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impose-command (continued) 

FATAL ERRORS: Xon-correctable syntax errors and invalid or ambiguous 
data. 
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let-command 

PURPOSE: To specify a macro. 

REFERENCE: Section £. 

FORMAT: let maovoname [ (p± s^ P2 ••• Sff-i Pn t s tf])] mean 
begi n 

d(piiP2i---iPN^ 
end 

where maovoname is an identifier, the p^ and sfe are arbitrary 
input elements other than dollar signs or parentheses, and 
d(pisp2i"-iPif) is an arbitrary sequence containing, in parti- 
cular, any number of occurrences of each p^, but not containing 
dollar signs or mismatched begin-end pairs (1 ^ fe ^ N, N ^ 1) . 

EXAMPLES: let ff mean begin bed /7fll.6*/ end 

let pr(f:x) mean begin print x using f end 
let dump(j from il to i2,) mean 

beg i n 

do(pr(ff:u(l,i,j),u(2,i,j)) for i=il step 1 until i2) 

end 

NORMAL OPERATION: The data base is modified so that during the process- 
ing of subsequent input the specified transformation will be performed. 
For transformation conventions, see Sections 9.1, 9.4, and 9.5. 

NORMAL RESPONSE: If maovoname was already active, a line of the form 

MACRO 'MACHONAME 1 NOW INACTIVE. 
Then, in any event, a line of the form 

INDICATED SUBSTITUTION FOR 'MACRONAME' WILL BE MADE. 

NON-FATAL ERRORS: Memory overflow. 

FATAL ERRORS: Non-correctable syntax errors. 
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1 i st-command 

PURPOSE: To cause listings of active identifiers to be printed. 

REFERENCE: Section 12.1 . 

FORMAT: list c±, c 2 , ■ ■ ■ , c N 

where each of the c fe is a class specifier, namely, one of the four 
words "curves", "variables", "macros", or "commands" (1 < k < N, 

N > 1). 

EXAMPLES: list macros, variables, curves 
list commands 

NORMAL OPERATION: The active identifiers, if any, for each c^, starting 
with c ± , are listed in the format shown below. Identifiers introduced 
by the user are listed in reverse chronological order. 

NORMAL RESPONSE: For each s^, either 

C k CURRENTLY ACTIVE: 'NAME ± , NAME 2 , ■■-, NAME M ' . 

or 

NO C k ACTIVE. 

depending on the current status of the indicated class. 

FATAL ERRORS: Syntax errors, e.g., misspelled class specifiers. 
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plot -command 



PURPOSE: To obtain or modify a graphical display. 



REFERENCES: Sections 1 and 7.2. 



FORMAT: plot 



Z 

aga i n 
ni 1 



[sans op -^ , op 2, 



°pn! 



off 

d i scard 
movie n 
tolerance n 

where z is an algebraic expression, n is a positive integer, and 
each of the opfo (I < k < N , N > 1) is an option designator, namely, 
one of the following: 



i 1 i nes 
(or j 1 i nes) 

spi kes 

base 

axes 
data 

d i fxyz 



zsca 1 e 
(or xyscal) 

xtrans 
(or ytrans) 
(or ztrans) 

s 1 ave 



causes lines connecting points of constant i (or j) 
to be omitted, except where needed to define surface 
or base-plane boundaries. 

causes vertical lines normally drawn at boundary dis- 
continuities from base plane to surface to be omitted 

causes outline of domain normally shown in base 
plane, plus spikes noted above, to be omitted 

causes indicator showing x- , y- , and z-directions 
(normally drawn at bottom left of screen) to be omit- 
ted 

causes numerical data indicating maximum and minimum 
values of x, y , and z (normally plotted at bottom 
right of screen) to be omitted 

causes same scale factor to be used for z -direction 
as is used for x- and zy-directions (normally s-coor- 
dinates are scaled independently, since their ranges 
may differ greatly from x,i/-ranges) 

causes new display to be constructed using same z- 
(or x,i/-) scale factor as was used in last display 

causes new display to be constructed using same x 
(or y~) (or z) translation constants as was used in 
last display 

causes displays to be plotted only on user's termi- 
nal, leaving second terminal free for another user 
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plot-command (continued) 

EXAMPLES: plot u(2, i ,j) 

plot atan (uy (l, i , j)/ux(l, i , j) ) sans slave 

plot again sans i 1 i nes , base 

plot off 

plot tolerance 8 

NORMAL OPERATION: If "plot z [sans ...]" is used, a command sequence of 
the form 

tal ly 

def i ne bed ' (z) ' = z 

is simulated and a data base established by evaluating "(z)" and lattice 
parameters "x" and "y" at all nodes. This data base is used to con- 
struct a (new) surface display satisfying the option list supplied by 
the phrase "sans op 1 , op 2 , ■■■, op N " , if the latter is included. Subse- 
quent displays of the same surface which satisfy differenct option lists 
can be generated from the same data base (hence at little additional 
computation expense) by using commands of the form "plot again [sans 
...]". Other commands have the following effects: 

(1) The command "plot nil" blanks the display tube and resets rota- 
tion angles, but leaves the user "signed on" to the display 
system. 

(2) The command "plot off" has the same effects as "plot nil", but 
signs the user off. 

(3) The command "plot cfiscard" causes the data base created as a re- 
sult of the last "plot ^''-command to be discarded, thereby free- 
ing memory for other purposes. 

(4) The command "plot movie n" causes the synchronized movie camera 
to expose n frames of film, both immediately and on each subse- 
quent use of push-button 4 (see Thornhill, et al [3], p. 2.4). 

(5) The command "plot tolerance n" causes a tolerance parameter used 
in constructing displays to be set ton. In general, increasing 
this parameter reduces picture quality, but increases the number 
of grid lines which can be displayed. For complicated surfaces 
based on lattices of many nodes, it is sometimes necessary to 
increase the tolerance parameter beyond its normal value of four 
in order to obtain a complete display. 
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plot -command (continued) 

NON-FATAL ERRORS: Illegal conditions recognized by the tally-command, 
memory overflow, references to undefined quantities, and various illegal 
conditions recognized by the display system. 

FATAL ERRORS: Non-correctable syntax errors and commands of the form 
"plot again [sans ...]" when no data base is active. 
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print-command 

PURPOSE: To produce printed output. 

REFERENCE: Section 7.1 . 

FORMAT: print e±, e2» • • • > &u [using bed X/X] 

where the e^ are algebraic expressions (1 5 k £ N, N ^ 1) , f is 
a MAD- language format specifier, and X is any single character 
that does not appear in /. 

EXAMPLES: print x, y, ux(l,-3,5), uy(2,7,9) 

print atan(y/x), 8*sxx-3.67 using bed /2f9.W 

NORMAL OPERATION: The e^ are evaluated and printed as indicated below. 

NORMAL RESPONSE: If the phrase "using bed X/X" is included, values are 
printed in accordance with format specifier /. Otherwise, they are 
printed in accordance with the format specifier "5(lpel6.7)*". See Sec- 
tion 7.1 for format restrictions and examples. 

NON- FATAL ERRORS: Memory overflow and references to undefined parame- 
ters . 

FATAL ERRORS: Non-correctable syntax errors, including illegal format 
specifiers. 
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quit-command 

PURPOSE: To return control to the CTSS supervisor. 
REFERENCES: Section 2, 11, and 13.2 . 
FORMAT: quit 

NORMAL OPERATION: Any disk file open as a result of a previous read- or 

record-command is closed and control is returned to CTSS via the A-core 

+ 
entry DORMNT. If control is passed back to EPS, e.g., by means of a 

CTSS start-command, processing of any input remaining in the command se- 
quence of which the quit-command was a part is continued in the normal 
manner. Thus by typing "quit a$", then a CTSS save-command, it is possi- 
ble to create a "saved" version of EPS which will process the input se- 
quence a immediately on being resumed. 

NORMAL RESPONSE: For any disk file that is closed, a message of the 
form 

FILE NAME-L NAME 2 HAS BEEN CLOSED. 

Then the usual CTSS time-usage message. 



+ See [2], sec. AG. 6. 01. 
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read -command 

PURPOSE: To initiate the reading of an input file or to change an input 
file's status from active to inactive or vice versa. 

REFERENCE: Section 11.1. 



FORMAT: read 



name i name 2 
swi tch 



where namey and name2 are identifiers, each with six or fewer 
characters . 

EXAMPLES: read alpha beta 
read swi tch 

NORMAL OPERATION: If the first form, i.e., "read namey name 2 " is used, 
any input file already open as a result of a previous read -command is 
closed, then name-y name 2 is opened and placed in active status. As a 
result, after the processing of the command sequence containing "read 
narney name 2" is completed, subsequent input is taken from namey name 2 ■ 
The reading of this file continues until 

(1) another command of the form "read name-y name2 u is encountered, 
in which case the procedure described above is repeated; 

(2] a quit-command is encountered, in which case the file is closed 
and control is returned to CTSS; 

(3) an end-of-file is encountered, in which case the file is closed 
and the user is requested to supply further input from the con- 
sole; 

(4) a "read switch" is encountered, in which case the file is placed 
in inactive status and, after processing of the command sequence 
containing "read switch" is completed, the user is requested to 
supply further input from the console; or 

(5) an error is detected, in which case the file is placed in inac- 
tive status and the user is immediately requested to supply fur- 
ther input from the console. 

If the reading of the file is interrupted because of condition (4) or 
(5), a "read switch" supplied from the console at any later time will 
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read -command (continued) 

cause reading to be continued from where it left off. 

NORMAL RESPONSE: For the command "read name-y name^ y the line(s) 

[FILE NAME' NAME' HAS BEEN CLOSED.] 
FILE NAME ± NAME 2 HAS BEEN OPENED. 

where the line in brackets is printed only if some input file name' 
name' was already open. For the command "read switch", either 

READ FILE HAS BEEN DEACTIVATED, 
or 

READ FILE HAS BEEN REACTIVATED. 

whichever is appropriate. 

NON- FATAL ERROR: Memory overflow. 

FATAL ERRORS: Invalid file names (including file names not in the 
user's directory), input files with improper format (see Section 11.1), 
and "read switch" commands when no input file is open. 
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record -command 

PURPOSE: To initiate the copying of console I/O into an output file or 
to change such a file's status from active to inactive or vice versa. 

REFERENCE: Section 11.2. 



FORMAT: record 



name-^ name 2 
swi tch 



where name-^ and name^ are identifiers, each with six or fewer 
characters . 

EXAMPLES: record gamma delta 
record swi tch 

NORMAL OPERATION: If the command "record namey name^' is used, any 
output file already open as a result of a previous record-command is 
closed, then name-^ name2 is opened and placed in active status. As a 
result, subsequent user input and program output are copied into the 
file name-^ name^, generating a transcript of all console I/O. (N.B., 
if name -1 name 2 refers to an existing file, information is appended to 
that file, not written over information already there.) The writing of 
a transcript continues until a quit-command or a "record switch" is 
given, the latter serving, in general, either to make an active output 
file inactive or to make an inactive file active. 

NORMAL RESPONSE: For the command "record name^ name 2" > the line(s) 

[FILE NAME' NAME^ HAS BEEN CLOSED.] 
FILE NAMh\ NAME 2 HAS BEEN OPENED. 

where the line in brackets is printed only if an output file name' name' 
was already open. For the command "record switch", either 

RECORD FILE HAS BEEN DEACTIVATED, 
or 

RECORD FILE HAS BEEN REACTIVATED. 

whichever is appropriate. 
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record -command (continued) 

NON- FATAL ERROR: Memory overflow. 

FATAL ERRORS: Invalid file names, track quota overflow, and "record 
switch" requests when no output file is open. 
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rel ax-command 

PURPOSE: To initiate or continue the successive relaxation of finite- 
difference equations. 

REFERENCES: Section 6, especially 6.2 , and Section 8. 

FORMAT: relax 

NORMAL OPERATION: If no difference equations exist, a form-command is 
simulated. Otherwise, it is assumed that the existing equations are to 
be maintained. Current values for relaxation factor "omega", tolerance 
parameter "delta", and maximum iteration index "limit" are used to up- 
date various internal parameters. Then relaxation is begun, using the 
current solution matrix as a starting point. During each iterative 
pass, corrections to the various components of the solution matrix are 
introduced in the same order as that in which difference equations are 
generated by the form-command. Thus all "n" dependent variables are 
treated at any given node before another node is considered, and, in 
terms of the lattice configuration in i ,j -space, nodes are processed 
from left-to-right, bottom-to-top. Iterations are continued until one 
of the two termination conditions associated with parameters "delta" and 
"limit" is satisfied, unless an arithmetic overflow is detected (gener- 
ally a result of numerical instability) or the user gives a single in- 
terrupt signal (see Section 13.1). 

NORMAL RESPONSE: A line of the form 

RELAXATION TERMINATED AFTER N PASSES. MAX SOLN CHANGE C. 

where N is the number of iterative passes completed and C is the maximum 
absolute change introduced in any component of the solution matrix dur- 
ing the /l/ th pass. 

NON-FATAL ERRORS: Illegal conditions detected as a result of simulated 
form-commands and arithmetic overflows. 
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remove -command 

PURPOSE: To delete macro definitions. 

REFERENCE: Section 9, especially 9.6 . 

FORMAT: remove begin maaronarne-^, macroname 2 , ... , macvoname^ end 
where the macvoname^ are identifiers (1 - k £ N, N > 1) . 

EXAMPLES: remove begin prescribe end 
remove begin pr, dump end 

NORMAL OPERATION: Each rnacronameu found actually to name a macro causes 
the deletion of that macro. 

NORMAL RESPONSE: For each macro actually deleted, a message of the form 
MACRO 'MACRONAME' NOW INACTIVE. 
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ret ry- command 

PURPOSE: To reinitiate processing of a command sequence saved as a 
result of a non-fatal error. 

REFERENCE: Section 14. 

FORMAT: retry 

NORMAL OPERATION: The most recently saved command sequence is taken 
from the command stack and retried. 

FATAL ERRORS: ret ry-cornmands given when no saved command sequences 
exists to be retried. 
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rev i ew-command 

PURPOSE: To determine current contour, parameter, or macro definitions, 
REFERENCE: Section 12.2. 



FORMAT: review 



curve 
var i abl e 
macro 



begin d 1 , d 2 , ■ • ■ , <% end 



where: (1) if the command "review curve ..." is used, each d^ is 
of the form 

ourvename [along i^ ,J]_, i 2 ,j 2 , ^3^3^ 
where ourvename is an identifier and the i^,J£ are algebraic ex- 
pressions (1 < I < 3); (2) if the command "review variable ..." 
is used, each dy_ is of the form 

paramname [ (s ]_ , s 2 , •■•> s^)] 
where paramname is an identifier and the s t are algebraic expres- 
sions (1 - £ - M > M - !); or ( 3 ) if the command "review macro ..." 
is used, each dy_ is of the form 

maoroname 
where macvoname is an identifier (I < k < N , N ^ I) . 

EXAMPLES: review curve begin zeta along 0,0, imax.O, imax,jmax end 
review variable begin x, y, a(2,2), h end 
review macro begin pr, dump end 

NORMAL OPERATION AND RESPONSE: In all cases, the dy_ are processed from 

left to right. 

If the command "review curve ..." is used, lists of grid coordi- 
nates and segment numbers are printed for each dy_. When a dy has the 
simple form ourvename , all nodes of contour ourvename are listed. When 
a dy has the qualified form "ourvename along i 1 ,j 1 , -i-2'^2' ^'^V' only 
the nodes of ourvename lying on the arc running from i 1 ,d 1 through i 2 ,j 2 
to £35,73 are listed. 
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rev i ew- command (continued) 

If the command "review variable ..." is used, a definition, value, 
or subscript range is printed for each dfc. For dfc which refer to unsub- 
scripted quantities or specific components of subscripted quantities, 
actual definitions (for quantities specified by def i ne-commands) or 
values (for quantities specified by set-commands) are printed. For dfa 
which refer to subscripted quantities, but do not include a full set of 
subscripts, current subscript ranges are printed. 

If the command "review macro ..." is used, the prototype call and 
corresponding definition for each dfc, i.e., each macroname , is printed. 

See Section 12.2 for examples. 

NON-FATAL ERRORS: References to undefined quantities among the grid 
coordinates i£,J£ of a "review curve ..." command or among the sub- 
scripts sp of a "review variable ...." command. 

FATAL ERRORS: Non-correctable syntax errors. 
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rotate-command 

PURPOSE: To cause surface displays to be pre-rotated by discrete amounts 

t 
before being transmitted to the display system. 

REFERENCE: Section 7.2. 

FORMAT: rotate e-^ about ai , e2 about a 2> •••> e N about a N 

where the efc are algebraic expressions and each afc is an axis 
specifier, namely, one of the symbols "x" , "y", or "z" (1 < k < N, 

N > 1) . 

EXAMPLES: rotate 180 about x 

rotate phi about z, 2»phi about y 

NORMAL OPERATION: The phrases "e^ about afc" are processed from left to 
right. Each efc is evaluated and interpreted as an angle (in degrees) 
about the afc-axis through which surface displays are to be pre-rotated. 
The rotation matrix associated with this operation is reset to the iden- 
tity matrix when a command such as "rotate about x" is received. 

NORMAL RESPONSE: The current display, if any, is replotted with the 
specified pre-rotations . Subsequent plot-commands produce new displays 
with the same pre-rotations until another rotate-command is received. 

NON-FATAL ERRORS: References among the e^ to undefined quantities and 
illegal conditions recognized by the display system. 

FATAL ERRORS: Non-correctable syntax errors. 



'Real-time rotations (controlled by the plexiglass globe) are, in 
effect, added to the pre-rotations introduced by the rotate-command. 
This command is not needed often, but in certain applications is useful 
for obtaining specific views. For example, by typing "plot nil rotate 
90 about x plot u(2,i,j)", the user will obtain immediately the view 
of surface "u (2 , i , j ) " seen by looking down the y-axis. Picture quality 
is generally superior when rotated views are obtained in this manner, 
for round-off errors introduced by the real-time hardware are avoided. 
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set-command 

PURPOSE: To specify algebraic parameters that are to be treated as 
constants . 

REFERENCES: Sections 3^5, 4.2 , 4.3 , 5.3 , 5.4, 6.2, 8.2 , 10.2, and 10.3. 

FORMAT: set Pu = p 12 = ••• = P^ = e l> P21 = P22 = ••• = 
P2M 2 = e 2> ■•■' Pffl = Pff2 = •'• = PNM N = e N 
where the py^ are parameter designators and the e^ are algebraic 
expressions (1 < k < N, N > 1; 1 < I < Mfc, % > 1) . 

EXAMPLE: set maxph i=pi= i t»atan (l) , z(imax)=pi/3 

NORMAL OPERATION: Argument groups containing more than one "=" are re- 
written as shown in the description of the def i ne-command, producing a 
revised command of the form 

set Pl = e x , p 2 = e 2 , ••-, Pl = e L 

For each assignment pfc = e^ of the revised command, the value of expres- 
sion efc is computed and the data base arranged so that in the future 
this value is fetched each time the value of p^ is required. Subscript- 
ed parameters are, again, treated as with the def i ne-command. Thus if 
the command of the example is given when "imax" is 32, the quantity 
"z(32) u is assigned the value ''pi/3". Note, since processing is from 
left to right, that "pi'' is known when the evaluation of "pi/3" is initi- 
ated. 

SIDE EFFECTS: Previous definitions are deleted and subscript ranges are 
modified as needed. 

NON-FATAL ERRORS: Memory overflow and references to undefined quanti- 
ties . 

FATAL ERRORS: Non-correctable syntax errors. 



■164- 



st op -command 

PURPOSE: Inserted before an if-phrase, to insure that any preceeding 
command is executed before the if-phrase is interpreted. 

REFERENCES: Sections 9^2, 10.2, and 10.3. 

FORMAT: stop 

NORMAL OPERATION: No action--control is returned to the EPS supervisor 
immediately. 
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d+ 



tal ly-comman 



PURPOSE: To close any bounding arcs not yet closed and to ready the 
portion of the data base which represents lattice structure in ^j'-space 
for the generation of finite-difference equations or for the construc- 
tion of graphical displays. 

REFERENCE: Section 3. 

FORMAT: tally 

NORMAL OPERATION: If a ta 1 1 y- command has been executed previously, 
either as a result of explicit user input or as a result of an earlier 
form- or plot-command, and if no append- or delete -command has been exe- 
cuted since that time, no action is taken. Otherwise, the following 
operations are performed. First, for each bounding arc cuvvename not 
yet closed, the command 

close ouT'venarrie 

is simulated. Then the node list for each contour is reordered, if nec- 
essary, so that the nodes are listed in clockwise (counterclockwise) or- 
der if the contour represents an external (internal) boundary. Finally, 
a check is made to insure that, in i ,j -space, no contour crosses over 
itself or over any other contour. At the same time, certain bookkeeping 
operations, including the determination of the total number of nodes in 
the lattice, are performed. 

NORMAL RESPONSE: None if no action is taken. Otherwise, a line report- 



Explicit use of the tally-command is unnecessary if input proce- 
dures described elsewhere are observed, for the command is simulated 
automatically whenever the operations it performs become essential to 
further processing, i.e., whenever a form- or plot-command is executed. 
Note, however, that it is possible to specify lattice structure in i,j- 
space by typing the relevant append-command (s) , then, in lieu of the 
usual close-command(s) , a single tally-command. This not only may save 
a small amount of typing, but also indicates immediately whether or not 
a lattice has acceptable topology. 
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tally-command (continued) 

ing the number of nodes N in the lattice as follows: 
POINT TALLY IS N. 

NON-FATAL ERRORS: Use of the command when no bounding contours are 
defined, when close-commands are needed but cannot be successfully exe- 
cuted, and when contours exist which cross over themselves or other con- 
tours . 
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16 . Summary 

To orient the beginning user and to provide the person who has 
used EPS in the past with a brief review, important input forms and 
notational conventions are outlined below. Items (3) through (9) con- 
stitute a check list which can be followed explicitly when solving line- 
ar boundary-value problems with no complicating factors such as free 
surfaces or undefined parameters. 

(1) To establish a link to EPS, type 

link norne^ saved ml4l6 cmfl04 

where nctme^ is "eps-p" if the ESL display is to be used for graph- 
ical output, or simply "eps" if otherwise. See Section 1. 

(2) To begin a run , type 

resume name^ 

and wait for the message "PROCEED:". In supplying subsequent 
input, remember that it is necessary to signal the completion of 
each unit of input explicitly by typing a "$" before giving a 
final carriage return. See Section 2. 

(3) To define lattice structure in i,j -space , type a sequence of the 
form 

append i±,3i, ^2»J2» •••» ^N'J'n to ^urvencane 
close ourvename 

for each bounding contour aurvename , where the ifctdk. 8* ve tne 
contour's vertices in the i ,j -plane in either clockwise or coun- 
terclockwise order. See Section 3, especially 3.3 and 3.4. 
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(4) To specify a_ mapping algorithm , type 

define x=f, y=g 

where f and g either explicitly or implicitly (by referencing 
auxiliary parameters) indicate how nodal values "x" and "y" for 
the independent variables depend on grid coordinates "i" and "j". 
If auxiliary parameters are introduced, specify functional depen- 
dencies by means of additional def i ne-commands , constants by means 
of appropriate set-commands. See Section 3, especially 3.5. 

(5) To specify governing partial differential equations , type set- 
and/or def i ne-commands designating parameters n, a{k,l), b(k,l), 
..., g(k,l), h(k) corresponding to n, a\j^, bfc£, ..., gfci, hfc in 

Yl 

1=1 Zx 9y (-^ 

+ f^T ( d kl ~^ + e kl^+ f\il u O + 9kl u l) = h k 

[k = 1,2,... ,n] 

Be certain that the fit-function is applied as shown in Section 
4.2 when equation parameters vary spatially in a continuous man- 
ner. 

(6) To number boundary segments corresponding to different boundary 
conditions, type commands of the form 

impose m along ii,j lt i2>3 2> ^3"? 3 

where m is a segment number to be associated with the arc running 
from node 'i\>3\ through ^2'3l to "^ 3 * J 3 • See Section 5.2. 

(7) To specify boundary conditions for a given segment m, type set- 
and/or def i ne-commands for parameters p(k,l,m), q{k,l ,m) , and 
r{k,m) corresponding to pfc^p Qklm' an ^ r km ^ n 
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n 

I ( Pkbn F l + VkJbn u l> " r fem W 

1=1 

[k - 1,2 n] 

where the F£ are first-order functions of the dependent variables 
U£ given in terms of the coefficients of governing system (1) as 
follows: 






(7) 



[fe. - 1,2,. ..,«] 

Minor restrictions on the Pfe^ and other details are given in Sec- 
tion 5.3. 

(8) To request the generation of finite-difference equations after 
data has been introduced as shown in items (3) through (7) , or 
modified as indicated in item (12), type 

form 

(9) To obtain a solution for the finite-difference equations , type a 
set-command giving values for relaxation factor "omega", tolerance 
parameter "delta", and maximum iteration index "limit". Then ini- 
tiate successive relaxation by typing 

relax 

If necessary, relaxation may be continued by means of additional 
re lax-commands. For information regarding the selection of values 
for "omega", etc., see Section 6.2. 

(10) To obtain printed output , type sequences such as 
pr i nt e^, e2» • • ■ > e N 
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do(print e^, «£» •••» e fl 

for \=ii step Si until iy) 

do(do(print e^, ei* -•>> e# 
for T«=^2. step 5t until i 2 ) 
for j=j'i step Sj until j' 2 ) 

where the e^ are algebraic expressions which may designate arbi- 
trary algebraic combinations of the solution properties u(k,i,j), 
ux.(k,i,j) , etc., as defined below in Table 2. See Section 7.1. 

(11) To obtain graphical output in the form of a surface display, type 

plot z 
where z is an algebraic expression. See Sections 7.2 and 15.2. 

(12) To modify lattice structure in i, ,j -space , type delete-, append-, 
and close-commands as explained in Section 8.1. To modify param- 
eter definitions , simply type more set- and/or def ine-commands. 
Similarly, to modify boundary segment numbers , type additional 

i mpose- commands . 

(13) To solve problems with nonlinearities or other features demanding 
iterative use of the system, combine items (8), (9), and (12), as 
explained in Section 10. 

(14) To avoid excessive typing , define appropriate macros by means of 
commands such as 

let maeroname{p± &\ Vl ••• S N-1 VB> me 311 
begin 

end 
where the p^ are formal parameters, the sfc are separators, and 
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d(p^;p2i • • • iViP i s a substitution to be made when calls of the 
indicated form are encountered. See Section 9. 

(15) To use disk input or to cause console I/O to be copied into a_ disk 
file , type read- and record-commands, as explained in Section 11. 

(16) To determine what identifiers are in use or to verify that defini- 
tions are correct , type list- and review-commands, as explained in 
Section 12. 

(17) To interrupt the machine in order to prevent or curtail the pro- 
cessing of erroneous input, press the interrupt button once. To 
terminate a run, type 

qu i t 

or press the interrupt button twice. See Section 13. 

Error messages resulting from improper input, etc., are generally self- 
explanatory, and in many cases indicate appropriate corrective action. 
However, if questions arise, consult Sections 14 and 15. 

The tables which follow summarize most of the important elements 
of the EPS input language. 
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Table 1. Summary of basic command forms 



f 



append i 1 ,j 1 , t 2 ,j' 2 > ■■•» *Wff [ to ^wename] 

close [aurvename] 

define Pll = p 12 = ... = p^ = e lt p 2 i = P22 
P2W 2 = e 2> •••' P/l/1 = PN2 = • 



= P/l% - e N 



expand begin p end 



form 



i gnore 



impose m 



along i 1 ,j 1 , [i 2 ^2^ ^3»J3 
at ix.J! 



[of curyename] 



1 e t maovoname [ (p ]_ s ^ p 2 
beg i n 

d(p±;p 2 ;- ■■;pu ) 

end 



s ^-l Ptf [ s /l/])] mean 



list c 1( c 2 . • • • . a N 

[sans op 1 , op 2 , 



plot 



2 

aga i n 
nil 



off 

d i scard 
mov i e n 
tolerance n 



opffl 



print e 1 , e 2 , ..., e^ [using bed A/A] 



f 



See sec. 15 for detailed descriptions, 
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Table 1 (continued) 



qu i t 



read 



name-^ name 2 
swi tch 



record 



name-^ name 2 
swi tch 



relax 



remove begin macvoname-^ , macroname2, •••, maoroname-^ end 



retry 



rev 1 ew 



curve 


var iable 


macro 

j 



begin d]_, ^2 > • • • » ■% enc ' 



rotate e^ about a-^, e2 about 0:2, ••■» <3/y about a^ 



set p-Q 

P2M 2 


= Pl2 = • 
= e 2 , •• 


• • ~ PlM ± ~ e l' 
• > PNI = PA/2 = 


P21 ~ P22 ~ • • 

• • • = PNM N = e N 


stop 








tal ly 
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Table 2. Summary of built-in functions 



Call 


Value 


Restrictions 


xs 


local boundary slope dx/d& 


use only in definitions of 
boundary condition parameters 


ys 


local boundary slope 9y/9-4 


p(k,l,m), q(k,l,m) and r{k,m) 
(see Section 5.3) 


octant 


lattice sector number 


use only in definitions of 
boundary condition parameters 
and governing equation param- 
eters a(k,l), b(k,l) , etc. 
(see Section 4.2) 


fit (a) 


interpolated or internodal 
value of z needed by dif- 
ference equation generator 
when z is to be treated as 
a continuous variable 


none- -produces the nodal val- 
ue of z for the current val- 
ues of i and j when differ- 
ence equations are not being 
generated (see Section 4.2) 


u(k,i,j) 


current approximation for 
UJ, at node (i,j) 


none (see Section 10 for pro- 
per use when solving nonline- 
ar equations) 


ux{k,i,,j) 


current approximation for 
dufc/dx at node (i,j) 




uy(k,i,j) 


current approximation for 
dufc/dy at node (t,j) 


use only in conjunction with 


uxx(k,i,j) 


current approximation for 
3 2 u/;./8x 2 at node (i,j) 


requests for output or with 
commands which cause problem 
data to be revised; do not 


uxy (k,i,j) 


current approximation for 
9 u^/9rr9i/ at node (i,j) 


use directly in definitions 
of boundary condition or gov- 
erning equation parameters 


uy x(k,i,j) 


current approximation for 
9 2 u^/9z/9x at node (i,j) 


(see Sections 7 and 10) 


uyy (k,i,j) 


current approximation for 
9 2 u^/9y at node (i,j) 




f]ux{k,i,j) 


current approximation for 
Ffc at node (i,j) 
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Table 3. Summary of operators 



t 



Oper- 
ator 


Type 


Prece- 
dence 


Usage 


Interpretation 




si n 






s i n a 


sine of a 






cos 






cos a 


cosine of a 


a in radians 




tan 






tan a 


tangent of a 






atari 






atan a 


arc tangent of a, result in radians 




log 


unary 


10 


log a 


natural logarithm of a [a > 0) 




exp 






exp a 


natural exponential of a, i.e., e 




sqrt 






sqrt a 


square root of a (a ^ 0) 




abs 






abs a 


absolute value of a 




+ 






+a 


plus a 




- 






-a 


minus a 




sign 


diadic 


9 


a sign b 


a with the sign of b 


power 


diadic 


8 


a power b 


a to the power b (a £ if b non- 
integral) 




/ 


diadic 


7 


a*b 
alb 


a multiplied by b 
a divided by b 




+ 


diadic 


6 


a+b 


a plus b 




- 






a~b 


a minus b 




grt 






a grt b 


1 if a > b, otherwise 




geq 






a geq b 


1 if a ^ b, otherwise 




les 
leq 


diadic 


5 


a les b 

a leq fc 


1 if a < b, otherwise 
1 if a ^ b, otherwise 




eq 1 






a eq] b 


1 if a = b, otherwise 




neq 






a neq b 


1 if a 4 b, otherwise 



















In expressions involving multiple operations, operators of great- 
er precedence are executed first unless parentheses indicate otherwise. 
Adjacent diadic operations of equal precedence are executed from left to 
right. Thus a/b*o is equivalent to {a/b)*s. Sequences of unary opera- 
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Table 3 (continued) 



Oper- 
ator 



not 



and 



or 



exor 



eqv 



Type 



unary 



diadic 



diadic 



diadic 



Prece- 
dence 



Usage 



not a 



a and b 



a or b 
a exor b 



a eqv b 



Interpretation 



1 if a = 0, otherwise 



1 if a ± and b + 0, otherwise 



1 if either a ± or b =f 0, or both, 

otherwise 

1 if either a t or b f 0, but not 
both, otherwise 



1 if a = and b = or if a ± and 
fo ^ 0, otherwise 



tors are, of course, executed from right to left, i.e., "exp sqrt abs a" 
can only mean "exp (sqrt (abs a))". 

Note, incidentally, that "exp", "sqrt", etc., are literally opera- 
tors in EPS, not built-in functions as in most compiler languages. Op- 
erands of such operators need not be enclosed in parentheses unless par- 
entheses are needed to indicate order of operation. Thus "sine of a" 
can either be written "sin a", as in mathematics, or "sin(a)", as in 
other languages; on the other hand, "sine of the quantity a+b" must be 
written "sin(a+fe)". 
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